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ABSTRACT 


Optimal  age  replacement  policies  are  designed  to  cut  down  system  failures  and 
minimize  maintenance  cost.  By  scheduling  planned  replacements,  a  system  is  replaced 
at  age  <t>*  or  at  failure,  whichever  comes  flrst,  and  the  cost  of  rq>lacement  before  failure 
(planned)  is  less  than  the  cost  after  failure  (unplanned).  In  this  thesis,  the  distribution  of 
lifetimes  is  a  known,  increasing  failure  rate  phase  type  distribution.  To  fmd  the  optimal 
age  of  replacement,  the  parameters  of  the  underlying  phase  type  distribution  must  be 
estimated. 

An  optimal  age  sequential  estimation  procedure  is  developed.  In  particular,  the 
phase  type  distributions  parameters  are  estimated  using  a  matching  moments  nonlinear 
programming  approach.  Since  there  are  many  parameters  associated  with  phase  type 
distributions  and  the  distributions  include  matrix  exponential  terms,  the  parameters  are 
in  general  difficult  to  estimate.  A  specific  case  where  the  phase  type  distribution  has 
initial  probability  vector  a =(1,0,0)  is  studied  for  different  sample  sizes  and  compared 
with  a  similar  nonparametric  procedure. 
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I.  INTRODUCTION 


Normally  a  system  or  component  is  replaced  when  it  fails  and  C,,  the  cost  of 
rq)lacement  before  failure  (planned  rq)lacement),  is  often  less  than  Cj,  the  cost  after 
failure  (unplanned  replacement).  The  problem  of  determining  when  to  rq)air  and  when 
to  replace  failing  systems  is  the  concern  of  management  resources.  Inefficient 
management  due  to  the  use  of  nonoptimal  maintenance  policies  can  lead  to  signiAcant 
system  maintenance  cost.  In  general,  optimal  maintenance  policies  are  designed  to  cut 
down  the  number  of  system  failures  and  minimize  maintenance  costs.  In  this  thesis,  we 
study  the  problem  of  estimating  a  type  of  optimal  maintenance  policy  (the  optimal  age 
replacement  policy)  when  the  underlying  system  life  distribution  is  phase-type.  In 
particular,  we  consider  an  adaptive  estimation  scheme  where  the  estimated  optimal 
replacement  age  is  updated  each  time  the  system  is  r^laced. 

Maintenance  is  defmed  to  be  all  activities  taken  to  keq>  the  system  in  serviceable 
condition  or  to  bring  it  back  to  serviceability.  There  are  two  types  of  maintenance, 
corrective  and  preventive  maintenance  [Ref.  5;  pp.  1-8].  Corrective  maintenance  is  used 
after  a  failure.  This  does  not  necessarily  mean  that  such  action  has  not  been  foreseen. 
Preventive  maintenance  aims  to  reduce  the  probability  of  failure  and  includes: 

-Planned  (scheduled)  maintenance  in  which  specified  components  are  replaced 
(usually  at  regular  intervals).  The  maintenance  time  is  usually  based  on  component 
lifetime  distributions. 
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-Unplanned  (condition-based)  maintenance  in  which  the  decision  to  replace  or  not 
to  replace  is  made  according  to  the  outcome  of  a  diagnostic  study. 

The  policies  that  are  designed  to  reduce  the  number  or  the  probability  of  system 
failure  and  maintenance  costs  are  called  maintenance  policies  [Ref.  21:  pp.  1-3].  We 
concentrate  on  age  replacement  policies  where  a  system  is  r^laced  at  age  T  or  failure, 
whichever  comes  first.  For  this  problem  to  make  sense,  Q  must  be  less  than  Cj 
system  must  age  with  time  i.e.  the  failure  rate  of  the  system  must  be  increasing  with  age. 
It  is  not  wise  to  replace  the  equipment  before  failure  when  the  system  has  a  constant  or 
decreasing  failure  rate  such  as  when  the  system  failures  occur  according  to  an  exponential 
distribution.  In  this  case,  rq}lacement  will  not  reduce  the  probability  of  system  failure 
occurring  in  the  next  instance  of  time.  When  using  an  age  rq)lacement  policy  the 
question  always  asked  is,  ”At  what  age  should  the  equipment  be  replaced”.  If  the 
replacement  occurs  too  frequently  the  maintenance  costs  will  be  high.  If  replacement  is 
too  infrequent,  the  system  will  fail  more  often  than  necessary  and  again,  the  maintenance 
cost  will  be  too  high.  Thus  it  is  desirable  to  Hnd  an  "optimal"  r^lacement  age.  Here, 
the  optimal  age  is  defined  as  the  age  which  yields  the  minimum  long  run  expected 
maintenance  costs  per  unit  time. 

We  will  assume  that  the  maintenance  action  returns  the  system  to  the  good-as-new 
condition,  thus  the  same  services  are  provided  as  before  rq)lacement.  To  accomplish  this 
we  assume  that  the  systems  used  for  rq}lacement  have  lifetimes  that  are  indq)endent  and 
identically  distributed  according  to  a  phase-type  distribution.  The  class  of  phase  type 
distributions  is  large  and  includes  Exponential  and  Gamma  distributions  along  with 
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convolutions  and  mixtures  of  these  distributions.  This  fact,  along  with  the  fact  that  phase 
type  distributions  have  a  physical  inteipretation  make  them  particularly  well  suited  for 
modeling  system  lifetimes. 

The  optimal  rq)lacement  age  depends  on  the  underlying  phase  type  distribution. 
This  is  in  general  unknown  and  must  be  estimated.  Estimation  for  phase  type 
distributions  based  on  iid  lifetimes  has  been  studied  by  Neuts  and  Meier  (1980). 
Estimation  of  the  optimal  age  of  replacement  for  phase  type  distributions  is  new. 
However  both  parametric  and  nonparametric  estimation  based  on  iid  observations  have 
been  considered  by  Arunkumar  (1972),  Ingram  and  Schaeffer  (1976),  Bergman  (1977), 
and  Barlow  (1978).  Here  we  use  a  sequential  approach  for  estimating  the  optimal 
replacement  age.  At  each  replacement  the  estimate  is  updated  and  the  new  system  is 
subject  to  an  age  replacement  policy  based  on  the  optimal  age  estimated  so  far.  This  type 
of  sequential  approach  has  been  studied  parametrically  by  Oclay  (1990)  and 
nonparametricly  by  Bather  (1977),  Frees  and  Ruppert  (1985),  Aras  and  Whitaker  (1991), 
Wu  (1990)  and  in  a  decision  theoretic  framework  by  Glazebrook,  Bailey  and  Whitaker 
(1991). 

Phase  type  distributions  arc  described  in  Chapter  U.  The  sequential  estimation 
procedure  to  estimate  the  optimal  age  of  replacement  is  given  in  Chapter  m,  and  in 
Chapter  IV  we  present  results  comparing  the  sequential  estimation  procedure  assuming 
an  underlying  phase  type  distribution  with  a  nonparametric  procedure.  Conclusion  and 
recommendation  are  given  in  Chapter  V. 
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n.  PHASE  TYPE  DISTRIBUTION 


A.  GENERAL 

A  phase  type  distribution  is  defined  as  the  distribution  of  the  time  until  absorption 
in  a  finite  state  Maikov  process.  To  determine  the  distribution  of  the  absorption  time, 
consider  an  m+1  state,  continuous  time  Markov  chain  whose  infinitesimal  generator  Q 
has  the  form 


where  T  is  a  nonsingular  m  x  m  matrix  with  negative  diagonal  elements  and  nonnegative 
off-diagonal  elements,  T®  is  vector  of  length  m  with  nonnegative  elements  and 
0=(0,0,...,0).  Moreover, 

Te  -H  T®  =  0' 
where  e  =  (1,1,...,!)'. 

By  the  construction  of  Q  the  state  m-i-1  is  absorbing  and  all  other  states  are 
transient.  The  necessary  and  sufficient  condition  for  this  is  that  the  inverse  T*  exists 
[Ref.  15:  pp.  446].  Let  the  initial  probability  vector  of  the  Maikov  chain  be  given  by  ( 
<^+1 )  oe  +  =  1  and  a  being  the  m  dimensional  initial  probability  vector 

of  transient  states  such  that  0  <  ae  ^  1.  LetXbethe  time  until  absorption  into  the 
(m+l)*^  state.  The  probability  distribution  F(x)  of  the  time  until  absorption  in  the  state 
m-l-1  corresponding  to  the  initial  probability  vector  (a,  a^n+i )  is  given  by 
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F(x)  =  1  -  aexp(Tx)e  ,  for  x  ^  0  , 


(2.1) 


where  exp(A)  is  the  matrix  exponential  of  the  matrix  A  defined  in  Ref.  8. 

The  probability  distribution  F(x)  on  [0,  oo)  is  said  to  be  of  phase  type  (PH- 
distribution)  and  the  pair  (a,  T)  is  called  a  representation  of  F(x).  If  >  0  the 
distribution  F(x)  has  a  jump  of  height  at  x  =  0. 

On  (0,  oo),  the  probability  mass  is  distributed  according  to  a  density  given  by 

f(x)  =  oeexp(Tx)T°  ,  for  x  >  0  ,  (2.2) 

The  noncentral  moments  /Xj  =  EpC']  of  F(x)  are  all  finite  and  given  by 

=  (-l)‘i!(ar‘e)  ,  for  i  ^  1  [Ref.  15:  pp.  446]  .  (2.3) 


B.  COST  FUNCTION 

Let  Xi,  Xj,...  be  a  sequence  of  independent  and  identically  distributed  (iid)  positive 
random  variables  with  distribution  function  F.  The  sequence  X,,  Xj,...  rqrresents  the 
sequence  of  system  lifetimes  that  would  be  observable  if  the  system  were  r^laced  at 
failure.  Let  C,,  Cj  (Cj  >  C,)  and  0  be  the  respective  cost  of  planned  replacement; 
unplaimed  replacement  and  the  replacement  age.  The  long  run  expected  cost  per  unit  time 
can  be  verified  [Ref.  1:  pp.  87]  to  be 


C{^) 


c,F(<|))  qF(4>) 

f*F(u)  du 
Jo 


(2.4) 


where  F(u)  =  1  -  F(u)  is  the  survival  function.  A  sufficient  condition  that  guarantees  a 
unique  and  finite  0*  that  minimizes  C(0),  is  that  F  have  strictly  increasing  failure  rate. 
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When  the  distribution  F  is  phase  type  with  representation  (a,  T),  the  long  run  expected 
cost  function  is 


^  Cg  [l-cexp(g\t>)a]  +  q[aexp(2\|))e] 

/*aexp  (Tti)  edu 
0 


(2.5) 


C;  -  aexp(Jj>)e[C;-q] 
[exp(2\{i)  -  J]  e 


(2.6) 


where  I  is  an  identity  matrix. 

The  phase  type  distribution  does  not  necessarily  have  increasing  failure  rate,  so  the 
cost  function  C(<^)  does  not  necessarily  have  a  unique  finite  minimum.  In  the  m=3  case 
even  when  the  initial  probability  of  the  transient  states  is  fixed  the  minimum  can  be 
unique  and  finite  or  can  occur  at  infinity  as  shown  in  Figure  1 ,  where  a  =  (1 ,  0,  0)  and 
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In  Figure  2,  the  same  phenomena  is  shown  when  F  has  the  same  nonsingular  3x3  matrix 
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but  the  initial  probabilities  are  different. 
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REPLACEMENT  AGE 


the  same  a  and  different  T  =  T,  or  Tj 
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m.  THE  SEQUENTIAL  ESTIMATION  PROCEDURE 

In  theory,  any  distribution  of  a  nonnegative  random  variable  can  be  approximated 
arbitrarily  closely  by  a  PH  distribution  [Ref.  11:  pp.  1-2].  This  means  that  the  PH 
distributions  are  dense  in  where  .^is  the  set  of  distributions  with  support  on  [0,oo), 
The  paper  of  Johnson  and  Taaffe  (1988)  shows  that  the  finite  mixtures  of  Erlang 
distributions  and  PH  distributions  with  a  Coxian  representation  are  both  dense  in  Since 
both  families  are  subsets  of  the  family  of  PH  distributions,  this  implies  that  PH 
distribution  are  dense  in  .^[Ref.  10:  pp.  1-8].  Thus,  the  PH  distribution  can  be  used  to 
approximate  any  unknown  distribution  of  lifetimes  (lifetime  is  always  greater  than  or 
equal  0).  Another  feature  that  makes  PH  distributions  a  desirable  choice  to  model  system 
lifetimes  is  that  they  may  have  a  physical  interpretation.  The  absorbing  state  represents 
system  failure  and  the  transient  states  may  represent  different  levels  of  a  system’s  ability 
to  function. 

A.  PARAMETER  ESTIMATION 

There  are  many  ways  to  estimate  the  parameters  of  a  distribution  including  the 
method  of  maximum  likelihood  (MLE),  and  the  method  of  moments.  The  PH 
distributions  have  special  properties  which  make  estimation  difficult.  First,  the  number 
of  parameters  is  not  fixed,  it  varies  with  the  number  of  transient  states,  e.g.  a  PH 
distribution  with  2  transient  states  has  6  {Munmeters,  a  PH  distribution  with  3  transient 
states  has  12  parameters  etc.  Also  the  probability  distribution  and  density  function  include 
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the  exponential  of  a  matrix.  This  makes  it  hard  to  woik  with  the  likelihood  function.  In 
this  thesis,  the  moment  matching  method  is  used  to  estimate  the  parameters  of  a  phase 
type  distribution  but  the  MLE  method  and  its  associated  problems  are  also  discussed  . 
In  addition,  the  number  of  transient  states,  m,  is  assumed  to  be  known. 

1.  The  Maximum  Likelihood  Method 

Let  Xj,  Xj,...  represent  the  sequence  of  system  lifetimes,  where  Xj,  Xj,... 
are  iid  PH  distributions  with  representation  (a,  T).  After  N  observations  the  data 
available  for  estimation  wiU  be  (Zj  6,),(Z2,52),...,(2i^,6N)  where  Zi=min  (Xi,$*.i).  X;  is 
the  i*  lifetime  and  is  the  optimal  age  replacement  estimated  after  i-1  observations  and 
5i=l  if  Xi^$*.i  5i=0  otherwise.  Assuming  system  lifetimes  have  a  PH  distribution 

the  likelihood  for  the  replacement  ages  Z,,...,2^  and  types  of  replacement  6,,  82 . 5n 

is 

jr 

L{a,T)  =  JJ  [eexp  (TZ_£)  r*]  *Maexp  e] 

i-1 


[  -oexp  ( TZj )  re]  [  aexp  ( TZ^ )  e]  ^ 


It  is  too  difficult  to  maximize  this  likelihood  directly  by  differentiating  L(a,T)  and 
solving  for  the  MLE’s.  There  are  several  reasons  for  this.  The  likelihood  and  its 
derivatives  include  the  exponential  of  a  matrix  form  that  cannot  be  written  in  closed  form 
and  the  parameters  are  subject  to  numerous  constraints.  In  addition,  there  are  many 
likelihood  equations  due  to  the  number  of  parameters  (3  transient  states  will  have  up  to 
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12  equations).  An  alternate  approach  is  to  find  approximate  MLE’s  using  nonlinear 
programming  algorithms  i.e.: 

jr 

MAXLia.T)  =  JJ  [-oexp  (rz_£)  Xe]  [aexp  (TZ^)  e] 

i«i 

S.T.  a  ^0 
ae  =  I 

tii  <  0  for  i  =  l,2,...,m 

tjj  ^0  for  i  j 

J^i^jtjj  ^  -tii  for  i  =  l,2,.,.,m 

det(T)  0 

where  m  =  the  number  of  transient  states  of  the  PH  distribution.  Even  when  m  is 
assumed  to  be  known,  there  is  still  the  problem  of  approximating  the  exponential  of  a 
matrix  and  the  optimization  software  to  take  care  of  this  problem  is  not  available.  Thus, 
it  is  very  difficult  to  use  this  approach  to  estimate  the  PH  distribution  parameters. 

2.  The  Moment  Matching  Method. 

The  PH  distribution  is  a  complicated  distribution,  there  are  many  parameters 
and  there  are  difficulties  with  computing  the  exponential  of  a  matrix.  One  property,  the 
existence  of  T’,  implies  that  all  noncentral  moments  of  a  PH  distribution  are  finite.  The 
k*  moment  is  given  by 

=  (-l)^!ar‘e.  (3.1) 

Moment  matching  is  a  common  method  for  estimating.  Using  moment  matching  by¬ 
passes  the  problem  of  evaluating  the  exponential  of  a  matrix  [Ref.  12:  pp.3]. 
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Let  m  be  the  number  of  transient  states  of  the  PH  distribution.  There  are  m(m+l) 
parameters  that  need  to  be  estimated,  thus  the  m(m+ 1)  equations  from  the  moments  that 
need  to  be  solved  are 

p.1  =  -aT'e 

Aj  =  2!ar'e 


=  (-ir""*^(m2+m)!aT‘e, 
where  is  the  k*  sample  moment. 

This  method  is  still  inappropriate  due  to  the  large  number  of  equations  and 
the  fact  that  when  the  dimension  of  matrix  T  is  big  and  the  moments  are  of  high  order, 
substantial  amount  of  error  is  introduced  when  trying  to  solve  these  equations.  Johnson 
and  Taaffe  have  shown  that  the  moment  matching  method  and  nonlinear  programming 
approaches  can  be  used  to  estimate  the  parameters  of  a  PH  distribution.  They  match  only 
the  second  and  third  standardized  moments  [Ref.  11:  pp.  3-11]. 

Let  /Xi  be  the  i***  central  moment.  The  second  standardized  moment,  the  coefficient 
of  variation  (C),  defmed  as 

^  _  standard  deviation 
mean 

can  be  written  as 
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c  = 


1^1 


(3.2) 


The  third  standardized  moment,  the  coefficient  of  skewness,  (7)  is  defined  as 


Y 


(3.3) 


And  the  t*  standardized  moment  is  for  t  =  3,4,.,.  . 

The  reason  for  matching  standardized  moments  is  that  the  standardized  moments  do  not 
change  with  scale  changes.  So,  the  shape  of  a  PH  distribution  with  representation  (a,T) 
will  not  change  if  it  is  multiplied  by  a  nonnegative  number. 

The  nonlinear  programming  formulation  for  matching  C  and  7  to  a 
continuous  PH  distribution  is  given  by 

MIN.  w,(C-C(0))2  +  Wj  (y-7(0))2 
S.T. 

a  ^  0 
ae  =  \ 

tji  <  0  for  i  =  1,2,. ..,m 

tjj  ^  0  for  i  j 

I tij  ^ -tii  for  i  =  l,2,...,m 

det(T)  0 

where 

C(0)  =  the  second  standardized  moment  of  the  cumrat  solution. 

7(0)  =  the  third  standardized  moment  of  the  current  solution. 
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m  =  the  number  of  transient  states  of  the  PH  distribution. 


ty  =  the  elements  of  row  i  and  column  j  of  nonsingular  matrix  T. 

Wi.Wj  are  positive  weights  chosen  to  guide  the  search  and  w,  ^  Swj. 

The  moments  are  matched  when  the  objective  function  value  is  zero.  The  solutions  which 
match  the  target  moments  are  called  "moment-matching  solutions"  [Ref.  11:  pp.  5]. 

Even  with  fixed  dimension,  this  procedure  may  have  many  solutions  and  these 
solutions  are  unpredictable.  There  are  some  practical  points  that  need  to  be  made 
associated  with  the  properties  of  PH  distribution  and  the  problem  of  using  nonlinear 
programming  to  get  the  moment  matching  solutions: 

-  Different  combinations  of  initial  solutions  and  algorithm  parameters  may  lead  to 
different  moment-matching  solutions.  Also,  some  combinations  of  initial  solution  and 
algorithm  parameters  may  not  lead  to  a  moment-matching  solution  and  the  nonlinear 
programming  algorithm  may  not  converge. 

-  For  a  given  dimension  we  do  not  know  how  many  solutions  exist  or  how  to  guide 
the  search  toward  a  preferable  solution. 

-  The  moment-matching  solutions  do  not  guarantee  that  the  estimated  cost  function 
will  have  a  finite  minimum.  When  choosing  among  several  solutions,  a  solution  which 
gave  a  cost  function  with  a  finite  minimum  was  always  chosen. 

-  The  values  of  w,  and  W2  can  be  modified  to  guide  the  search;  generally  w,  ^  1, 
W2^  1  and  w,  ^  3w2.  The  magnitude  of  w,  and  W2  can  be  increased  to  get  more  or 
sufficient  accuracy  [Ref.  11:  pp.  6]. 
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-  If  the  feasible  solutions  are  known,  parameters  can  be  bounded  to  guide  the 
search.  However,  some  bounds  on  the  parameters  may  not  lead  to  a  moment-matching 
solution  or  to  a  solution  whose  cost  function  does  not  have  a  unique  finite  minimum 

-  User  interaction  is  often  necessary  in  obtaining  the  appropriate  or  preferable 
solutions. 

B.  THE  SEQUENTIAL  ESTIMATION  PROCEDURE 

Let  X,,X2, . ,Xn  be  the  sequence  of  lifetimes  of  the  system  and  {$*n}  be  the 

sequence  of  estimators  of  0*  where  is  the  estimator  after  N  replacements.  Then  after 
N  observations,  we  have  the  right  censored  data  (Zj,6,),(Z2,62),.-..,(Zn.5n)- 
where  Zj  =  min  (Xi,$Vi) 

Sj  =  1  if  Xj  :S  $*i.i  otherwise  5;  =  0  . 

Because  the  data  are  right  censored  the  usual  method  of  moments  approach  for  estimation 
needs  to  be  modified.  Rather  than  use  sample  moments  calculated  from  an  empirical 
distribution,  we  use  moments  from  a  nonparametric  estimator  of  F  based  on  the  censored 
data.  From  the  right  censored  data  after  N  observations,  the  procedure  to  compute  the 
estimators  {$’ }  is  developed  as  follows. 

1.  Use  the  right  censored  data  sequence  (Z„6j),(Z2,5i),  ...,(Zn,6n)  to  estimate 
moments  with  F  =  l-F  by  the  product-limit  estimator: 
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(3.4) 


F(t) 


n 


I  ^  ^  )  *(1> 

N-i+1 


where  Z(,),  Z^  are  the  order  statistics  of  Z,,  Z^,...,  Z^  and  6(,),  fi®,...,  6(N)  are 

ordered  according  to  the  ordering  of  Z^^,  Zp),...,  Z^^).  Then  calculate  probabilities  Pj, 
associated  with  Z^^,  by 

Pi  =  F(Z,0-f(Z(i)).  (3.5) 

Note,  since  F  has  discontinuities  only  at  Z^  where  =  1 ,  that  Pj  =  0  when  6(j)  =  0. 
Then  estimates  of  the  k*  moment  are  given  by 


The  estimate  of  the  2*^  standardized  moment  is 


r  -  -(gUl)^ 

m] 


and  the  estimate  of  the  3'**  standardized  moment  is 


(3.6) 


(3.7) 


{y/E[X^]  -iE[X])^)^ 

2.  Estimate  parameters  of  PH  distribution  (a,T),  by  matching  the  second  and  third 
standardized  moments  using  the  nonlinear  programming  approach. 

When  executing  the  nonlinear  programming  algorithm  using  the  model  described 
in  the  previous  section,  the  det(T)  0  constraint  is  Tcphced  by  a  constraint  on  the 
expected  lifetime.  The  reason  for  doing  this  is  that  we  could  not  formulate  an  algebraic 
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constraint  equivalent  to  the  constraint  det(T)  ^  0  .  The  problem  can  be  solved  by  a 
constraint  on  the  expected  lifetime,  i.e.  by  taking 

aT-^e  =  i![X]  <3  *9) 

where  E[X]  is  calculated  in  step  1.  This  constraint  will  make  the  search  for  preferable 
solutions  easier.  The  initial  solution  is  chosen  by  the  user.  Some  initial  solutions  may  or 
may  not  lead  to  a  moment  matching  solution.  Typically,  the  initial  vector  a  is 
(1/m, 1/m,. ..,l/m)or(l,0,...,0) and  an  initial  matrix  T  consists  of  tj  j  =  -1  andtjj  =  1/m 
for  i  j,  i  =  1,2,. ..,m.  In  this  thesis  the  initial  vector  a  is  (1,0,. ..,0)  and  an  initial 
matrix  T  is  taken  to  have  tjj  =  -1  and  tjj  =  1/m  [Ref.  11:  pp.  9]  . 

3.  Using  the  parameters  estimated  in  step  2,  the  estimated  cost  function  is 

d(^)  =  ~  &exp_(a^)a[q-Ci]  (3.10) 

af'Mexptf^)  -  I]o 

The  new  estimate  of  the  optimal  age  of  r^lacement  is  taken  to  be  the  0  that 
minimizes  (3. 10)  by  enumerative  search. 

4.  Compare  the  with  Xn+,  then  repeat  Step  1. 

A 

The  initial  solutions  for  matrix  T  and  vector  a  in  step  2  are  taken  to  be  the  previously 
estimated  values  of  T  and  a.  In  case  the  moment  matching  solution  cannot  be  met,  the 
initial  solution  in  step  2  will  be  taken  to  be  the  original  initial  solution  and  the  previous 
optimal  age  rq)lacement  is  used  for  stq)  3. 

The  rq)lacement  cost  for  i"*  system  is  C2  if  X;  <  otherwise  the  rq7lacement 
cost  is  C|.  The  actual  total  rqtlacement  cost  for  the  Hrst  N  systems  that  are  observed  is 
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(3.11) 


N 


and  the  total  operating  time  for  the  N  systems  is 

.  (3.12) 

1-1  i-l 

So,  the  actual  average  cost  (AAC)  after  N  r^lacement  is  computed  by 


AAC^  =  (3.13) 

In  this  thesis,  the  moment  matching  nonlinear  programming  approach  for  estimating 
parameters  of  PH  distributions  uses  the  GAMS  program  as  shown  in  .^rpendix  A  (for 
a  PH  distribution  with  3  transient  states).  The  Fortran  programs  are  written  to  estimate 
the  optimal  age  of  replacement,  second  and  third  standardized  moments  and  for  data 
preparation  as  shown  in  Appendix  B. 
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rV.COMPARlSON  WITH  NONPARAMETRIC  PROCEDURE 


A.  NONPARAMETRIC  PROCEDURE 

An  alternate  approach  for  estimating  is  to  estimate  nonparametrically ,  as  in  Aras 
and  Whitaker  (1991).  The  procedure  is  much  simpler  computationally  than  the  parametric 
procedure  described  in  the  previous  section  for  PH  distributions.  In  the  nonparametric 
procedure  after  N  replacements  the  product  limit  estimator  F  of  F  given  in  (3.4)  is  used 
to  estimate  the  cost  function  C(0)  as 


Cs(^)  - 


C2#((t>)  +  qF((|>) 

f  riu)  du 
J  0 


(4.1) 


^  — 

where  F  =  1  -  F  is  the  estimator  of  the  cdf  F,  the  estimator  of  then  found  by 
minimizing  In  general  one  would  expect  parametric  procedures  to  do  much  better 
than  nonparametric  procedures.  However,  with  the  numerical  difficulties  involved  with 
estimating  parameters  for  PH  distribution  and  the  fact  that  the  family  of  PH  distributions 

4, 

is  so  large  it  is  not  obvious  which  procedure  is  best.  The  criterion  for  comparison  is  the 
actual  average  cost  per  unit  time,  AAC^,  after  N  replacements. 
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B.  COMPARISON  OF  NONPARAMETRIC  AND  PARAMETRIC  PROCEDURE 


FOR  THE  PH  DISTRIBUTION 

Simulation  was  used  to  compare  sequential  estimation  based  on  PH  distributions 
with  nonparametric  sequential  estimation,  C,  =  100,  Cj  =  500.  System  lifetimes  were 
generated  by  an  increasing  failure  rate  PH  distribution  with  representation  (T4,,a4,) 
where, 


‘41 


-5  4  1 

I- 6  5 

II- 7 


a4i  =  (1.0, 0.0,  0.0) 


For  this  representation  the  long  run  expected  cost  per  unit  time  is  given  in  Figure  3. 

To  give  the  parametric  procedure  a  chanc'*  ’  c  first  25  system  lifetimes  are 
uncensored.  After  the  first  25  observations,  the  parameters  of  the  PH  distribution  are 
estimated  sequentially  as  in  Chapter  El  Section  B.  The  actual  average  costs  (AAC)  are 
calculated  and  compared  for  small,  medium  and  large  sample  sizes  with  planned  cost  C, 
=  100  and  unplanned  cost  Cj  =  500.  The  20  initial  data  sets  are  simulated  as  shown  in 
the  Appendix  A. 

The  programs  are  separated  into  5  parts.  The  first  one  is  written  in  GAMS  as 
shown  in  Appendix  A,  to  estimate  the  parameters  of  a  PH  distribution  by  the  nonlinear 
programming  approach.  The  rest  are  written  in  Fortran  as  shown  in  Appendix  B,  to 
estimate  the  optimal  age  of  replacement  based  on  the  estimated  PH  distribution  from  the 
GAMS  program;  to  simulate  the  system  lifetime;  to  prepare  the  right  censored  data;  to 
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estimate  the  expected  lifetime,  second  and  third  standardized  moments;  and  finally  to 
tabulate  the  results  and  to  prepare  the  initial  data  for  the  next  estimation. 


LONG  RUN  EXPECTED  COST  PER  UNIT  TIME  OF  PH  DISTN. 
WITH  representation  T41.  o41  AND  C1  =  100.  C2=500 


REPLACEMENT  AGE 


Figure  3  Long  run  expected  cost  per  unit  time  curve  of  PH  distribution  with 
representation  T^,  and  a4i 

In  this  simulation,  the  number  of  transient  states  is  fixed  at  3  and  the  initial 
probabUity  vector  a  =  (1,  0,  0).  By  fixing  a  =  (1,  0,  0)  we  reduce  the  number  of 
parameters  to  be  estimated  and  use  a  model  in  which  all  systems  start  in  the  same  state, 
when  states  are  thought  of  as  the  level  of  system  repair,  this  choice  better  reflects  the 
idea  that  replacement  systems  are  as  good  as  new.  Bounded  parameters  and  user 


21 


interaction  are  used  in  guiding  each  estimation  during  the  simulation  to  a  preferable 
solution. 

Tables  1,  2  and  3  summarize  the  simulation  results  of  the  actual  average  cost 
resulting  from  sequential  estimation  of  PH  and  the  nonparametric  procedure.  Also  given 
are  the  signed  ranks  of  the  differences  in  the  actual  average  costs.  These  are  used  in  the 
Wilcoxon  Signed-Rank  test  to  test  the  hypothesis: 

Ho:  E[AACph]  =  E[AACnonp] 

H,:  E[AACph]  <  E[AAC)40Np] 

by  using  T^  the  sum  of  the  positive  ranks  as  the  test  statistic  . 

Tables  1,  2  and  3  give  the  statistic  T^  as  35,  10  and  2  respectively.  When  compared  to 
the  value  in  Table  9  [Ref.  13:  pp.  780]  at  N=20  all  p-values  are  less  than  0.005  which 
implies  that  the  sequential  estimation  of  a  PH  distribution  is  better  than  sequential 
estimation  of  nonparametric  for  all  sample  sizes  (small,  medium  and  large).  The  bc\ 
plots  in  Figures  5  and  6  of  actual  average  cost  versus  the  number  of  replacements  at 
different  sample  sizes  show  that  the  actual  average  costs  decrease  as  the  number  of 
replacements  increase  and  that  AACpH  for  the  parametric  procedure  decrease  more  than 
the  AACnohp  foi*  nonparametric  procedure. 

A  second,  PH  distribution  was  chosen  to  compare  the  parametric  PH  procedure 
with  the  nonparametric  procedure.  The  second  PH  distribution  has  an  average  long  run 
expected  cost  function  that  is  more  shallow  than  the  first  PH  distribution.  It  has 
representation  T42  and  042 


22 


where 


-5  4  1 

=  1  -6  5  ,  ««  =  (1.0,  0.0,  0.0)  , 

12-7 


and  the  long  run  expected  cost  per  unit  time  is  given  in  Figure  4. 


Figure  4  Long  run  expected  cost  per  unit  time  curve  of  PH  distribution  with 
representation  T4}  and  042 

Tables  4,  5  and  6  summarize  the  result  and  give  the  statistic  of  93,  81  and  47 
respectively.  The  p- value  taken  from  table  9  [Ref.  13:  pp.  780]  indicated  that  the 
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parametric  procedure  is  not  better  than  the  nonparametric  procedure  for  small  and 
moderate  sample  sizes.  In  fact,  for  some  cases  (e.g.  runs  16  and  20  in  Table  4  and  runs 
1,  5  and  16  in  Table  5)  the  actual  average  cost  resulting  from  using  the  nonparametric 
procedure  can  be  quite  a  bit  less  than  the  actual  average  cost  resulting  from  the 
parametric  procedure.  This  is  due  to  the  fact  that  C(4>)  is  shallow  at  thus  even  though 
C(<f>)  may  be  a  reasonable  estimator  of  C(0),  the  variance  of  its  minimizing  value  is 
higher  than  for  the  previous  PH  distribution.  This  means  that  is  often  underestimated. 
Because  C(<t>)  increases  so  rapidly  to  the  left  of  </>*,  underestimating  can  increase  the 
actual  average  costs  dramatically.  Another  difficulty  is  that  the  estimated  Ci<l>)  may  be 
relatively  flat  around  This  causes  numerical  difficulties  with  the  nonlinear 
programming  algorithm.  For  large  sample  sizes  the  parametric  procedure  does  do  better 
than  the  nonparametric  procedure.  Here  C(^)  estimates  C(^)  more  closely  and  it  is  easier 
for  the  nonlinear  programming  algorithm  to  identify  the  correct  solution. 

In  Figures  5  and  6,  the  improvement  in  actual  average  cost  with  sample  size  is 
shown  for  both  the  parametric  and  nonparametric  procedure  for  the  PH  distribution  with 
representation  T4,  and  a4,.  As  expected,  both  procedures  improve  (actual  average  costs 
decrease)  with  sample  size.  However,  the  actual  average  cost  for  the  parametric 
procedure  decreases  faster  than  for  the  nonparametric  procedure.  Both  procedures  exhibit 
considerable  variability.  The  solutions  have  got  still  mixed  that  we  should  do  more  study 
on  another  representation  until  the  more  precise  solutions  come  up. 
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Table  1  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  SMALL  SAMPLE 


SIZES  (N  =  50)  OF  PH  DISTRIBUTION  WITH  T4,  AND  «« 


Run 

AACpH 

AACnon? 

AACkp 

DIFFERENCE 

RANK 

1 

690.095 

754.555 

759.142 

-64.46 

15 

2 

714.728 

775.712 

780.872 

-60.444 

13 

3 

700.342 

684.379 

688.578 

15.963 

4  * 

a 

666.014 

725.331 

738.987 

-59.317 

12 

5 

714.330 

757.528 

735.058  $ 

-43.198 

8 

6 

544.322 

621.505 

621.507 

-77.183 

19 

B 

721.505 

810.149 

810.149 

-88.644 

20 

8 

565.854 

617.458 

636.556 

-51.604 

11 

9 

676.511 

725.394 

702.676  $ 

-48.883 

10 

10 

698.572 

766.182 

809.629 

-67.61 

17 

11 

814.989 

767.586 

801.150 

47.403 

9 

12 

644.921 

667.553 

672.116 

-22.632 

5 

13 

761.516 

686.794 

679.892  ! 

74.722 

18  * 

14 

608.477 

601.215 

617.711 

7.262 

1  • 

15 

706.591 

767.782 

782.198 

-61.191 

14 

16 

629.294 

667.805 

682.719 

-38.511 

6 

17 

688.572 

755.475 

748.353  $ 

-66.903 

16 

18 

636.963 

624.132 

665.344 

12.831 

3  * 

19 

583.226 

625.161 

609.036  $ 

-41.935 

7 

20 

644.065 

654.691 

676.005 

-10.626 

2 

AACkf  =  Actual  average  cost  rq)lacing  when  failure  occurs 


DIFFERENCE  =  AACpH  -  AACnqnp 
♦  positive  rank 
T^  =  35  ; 
p-value  <  0.005 

$  AACrp  is  less  than  only  AACnonp 
!  AACu;  is  less  than  AACpH  and  AACnon? 
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Table  2  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  MEDIUM 


SAMPLE  SIZES  (N  =  100)  OF  PH  DISTRIBUTION  WITH  T4,  AND 


Run 

AACpH 

AACnonp 

AACrf 

DIFFERENCE 

RANK 

1 

648.127 

723.170 

732.263 

-75.043 

18 

2 

690.971 

786.554 

789.767 

-95.583 

20 

3 

649.886 

644.710 

683.933 

5.176 

2  * 

n 

696.858 

747.640 

755.982 

-50.782 

11 

5 

690.195 

688.749 

682.840  ! 

1.446 

1  * 

6 

577.742 

607.414 

610.850 

-29.672 

6 

B 

686.414 

766.365 

762.245  $ 

-79.951 

19 

8 

550.106 

612.046 

638.286 

-61.94 

14 

9 

637.268 

696.249 

689.552  $ 

-58.981 

12 

10 

684.253 

730.138 

786.920 

-45.885 

10 

11 

802.039 

790.698 

807.530 

11.341 

4  • 

12 

615.087 

687.719 

695.174 

-72.632 

17 

13 

677.302 

666.119 

674.830 

11.183 

3 

14 

586.838 

628.687 

638.451 

-41.849 

8 

15 

726.213 

749.507 

737.792  $ 

-23.294 

5 

16 

626.605 

694.232 

667.462  $ 

-67.627 

16 

17 

671.908 

732.053 

758.912 

-60.145 

13 

18 

632.954 

667.156 

727.346 

-34.202 

7 

19 

635.989 

679.399 

667.141  $ 

-43.41 

9 

20 

590.881 

657.064 

701.735 

-66.183 

15 

> 

average  cost  re 

placing  when 

failure  occurs 

DIFFERENCE  =  AACpH  -  AACnonp 
*  positive  rank 
T^  =  10  ; 
p-value  <  0.005 

$  AACrf  is  less  than  only  AACnon? 

!  AACitp  is  less  than  AACpH  ^^nonp 
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Table  3  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  LARGE  SAMPLE 


SIZES  (N  =  200)  OF  PH  DISTRIBUTION  WITH  T4,  AND  «« 


Run 

AACpH 

AACnonp 

AACpp 

DIFFERENCE 

RANK  1 

1 

622.560 

655.891 

655.865  $ 

-33.331 

7 

2 

658.640 

711.831 

708.966  $ 

-53.191 

11 

3 

623.823 

635.065 

671.611 

-11.242 

1 

a 

678.141 

749.866 

754.490 

-71.725 

16 

5 

674.548 

714.291 

707.926  $ 

-39.743 

9 

6 

585.435 

645.920 

647.927 

-60.485 

13 

a 

647.302 

729.159 

733.860 

-81.857 

20 

8 

581.711 

607.068 

649.386 

-25.357 

5 

9 

638.465 

678.078 

719.904 

-81.439 

19 

10 

655.709 

717.569 

766.020 

-61.86 

14 

11 

760.478 

779.362 

762.986  $ 

-18.884 

3 

12 

589.088 

649.545 

681.004 

-60.457 

12 

13 

696.081 

679.816 

712.761 

16.265 

2  • 

14 

602.032 

643.703 

649.832 

-47.8 

10 

15 

685.977 

767.133 

755.331  $ 

-81.156 

18 

16 

646.752 

723.719 

663.489  $ 

-76.967 

17 

17 

680.457 

713.146 

754.616 

-32.689 

6 

18 

674.300 

697.304 

740.394 

-23.004 

4 

19 

636.649 

674.845 

670.146$ 

-38.196 

8 

20 

588.468 

654.994 

670.194 

-66.526 

15 

■H||V 

VACbp  =  Actual  average  cost  replacing  when  failure  occurs 

DIFFERENCE  =  AACpH  -  AAC^onp 
*  positive  rank 
T*  =  2  ; 
p- value  <  0.005 

$  AACkf  is  less  than  only  AAC^^onp 
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Table  4  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  SMALL  SAMPLE 


SIZES  (N  =  50)  OF  PH  DISTRIBUTION  WITH  T42  AND  042 


Run 

AACpH 

AACnonp 

AACpp 

DIFFERENCE 

RANK 

1 

618.463 

561.050 

545.871  ! 

57.413 

16  * 

2 

604.098 

684.564 

686.650 

-80.466 

19 

3 

674.353 

664.878 

637.614  ! 

9.475 

6  •  1 

m 

605.528 

559.966 

575.315 

45.562 

15  *  1 

5 

580.054 

521.136 

535.683 

58.918 

17  *  1 

6 

547.534 

553.458 

556.857 

-5.924 

4 

B 

632.508 

625.711 

625.711 

6.797 

5  * 

8 

632.510 

647.089 

651.264 

-14.579 

7 

9 

490.099 

578.531 

566.695  $ 

-88.432 

20 

10 

629.150 

668.268 

670.786 

-39.118 

13 

11 

535.032 

536.771 

535.194 

-1.739 

3 

12 

605.358 

644.458 

645.689 

-39.100 

12 

13 

519.537 

552.440 

552.443 

-32.903 

10 

14 

532.469 

553.664 

569.868 

-21.195 

8 

15 

671.880 

672.507 

672.160$ 

-0.627 

1 

16 

645.125 

571.650 

576.024 

73.475 

18  * 

17 

495.932 

533.334 

539.149 

-37.402 

11 

18 

648.741 

677.909 

682.898 

-29.168 

9 

19 

577.580 

576.101 

600.152 

1.479 

2  * 

20 

560.629 

519.418 

532.841 

41.211 

14  *  1 

AACrf  =  Actual  average  cost  rq)laciiig  when  failure  occurs 
DIFFERENCE  =  AACpH  -  AACnonp 


*  positive  rank 
T+  =  93  ; 
p-value  >  0.05 

$  AACpp  is  less  than  only  AACnqnp 
!  AACrf  is  less  than  AACpH  and  AACnonp 
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Table  5  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  MEDIUM 
SAMPLE  SIZES  (N  =  100)  OF  PH  DISTRIBUTION  WITH  T42  AND 


Run 

AACpH 

^^ACnoNP 

AACrf 

DIFFERENCE 

RANK  1 

1 

573.310 

548.751 

539.868  ! 

24.559 

13  » 

2 

588.501 

633.798 

623.661  $ 

-45.297 

16 

3 

618.891 

626.148 

608.227  ! 

-7.257 

2 

a 

621.602 

590.532 

601.378 

31.070 

14  * 

5 

557.189 

521.040 

535.590 

36.149 

15  * 

6 

549.717 

562.877 

569.251 

-13.160 

5 

D 

651.656 

629.631 

630.823 

22.025 

10  * 

8 

615.798 

624.193 

628.347 

-8.395 

4 

9 

475.757 

574.707 

567.688  $ 

-98.95 

20 

10 

602.762 

626.458 

628.035 

-23.696 

12 

11 

535.426 

537.073 

549.888 

-1.647 

1 

12 

580.480 

598.790 

590.995  $ 

-18.310 

7 

13 

525.770 

597.208 

597.210 

-71.438 

19 

14 

518.669 

540.978 

553.372 

-22.309 

11 

15 

634.978 

627.571 

633.019 

7.407 

3  * 

16 

675.027 

616.239 

618.041 

58.788 

18  * 

17 

535.911 

585.289 

587.998 

-49.378 

17 

18 

617.194 

636.912 

645.975 

-19.718 

9 

19 

550.468 

564.214 

579.831 

-13.746 

6 

20 

1  539.256 

520.034 

531.189 

19.222 

8  * 

AACrp  =  Actual  average  cost  rqjlacing  when  failure  occurs 
DIFFERENCE  =  AACpH  -  AACnonp 
*  positive  rank 


T^  =  81 


9 


p-value  >  0.05 

$  AACkf  is  less  than  only  AAC,^onp 
!  AACpp  IS  less  than  AACpn  and  AACi^onp 
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Table  6  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  LARGE  SAMPLE 
SIZES  (N  =  200)  OF  PH  DISTRIBUTION  WITH  T42  AND  a42 


AACpH 

AACnonp 

AACkp 

DIFFERENCE 

556.373 

550.131 

544.302  ! 

6.242 

574.206 

616.553 

621.391 

-42.347 

578.172 

591.186 

606.516 

-13.014 

594.573 

579.915 

588.547 

14.658 

580.261 

551.976 

560.463 

28.285 

526.351 

558.604 

573.294 

-32.253 

632.238 

622.079 

624.275 

10.159 

612.338 

621.511 

610.932  ! 

-9.173 

460.437 

567.069 

563.305  $ 

-106.632 

557.424 

599.181 

602.122 

-41.757 

555.7’ i 

539.303 

551.773 

16.408 

56G.534 

584.128 

589.078 

-23.594 

576.939 

621.856 

621.857 

-44.917 

568.925 

595.619 

599.305 

-26.694 

562.097 

583.057 

595.668 

-20.960 

584.778 

604.211 

606.979 

-19.433 

549.979 

607.061 

610.878 

-57.082 

631.719 

628.295 

639.159 

3.424 

578.172 

577.936 

607.344 

0.236 

538.683 

530.002 

570.516 

8.681 

RANK 


3 


17 


kf  =  Actual  average  cost  r^lacmg  when 
DIFFERENCE  =  AACpH  -  AACnonp 
*  positive  rank 
T*  =  47  ; 

0.01  <  p-value  <  0.025 
$  AACkf  is  less  than  only  AAC^onp 
!  AACup  IS  less  than  AACp^  and  A^^Cj^qi^ 


ure  occurs 


500 


ACTUAL  AVERAGE  COST 

600  700  600 


Figure  6  The  box  plot  of  actual  average  cost  for  different  sample  sizes  of 
nonparametric  sequential  procedure  for  T4,  and  a4, 
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UENTIAL  estimation  PROCEDURE 


V.CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  thesis,  the  age  r^lacement  policy  has  been  considr.ied.  A  system  is  replaced 
at  the  time  of  failure  or  at  a  scheduled  replacement  age  whichever  comes  first.  The 
replacement  system  is  assumed  to  be  as  good  as  new.  The  objective  is  to  achieve  the 
minimum  long  run  expected  maintenance  cost  per  unit  time.  The  system  lifetimes  are 
assumed  to  have  PH  distributions.  Estimating  the  parameters  of  a  PH  distribution  is 
difficult  because  it  may  have  many  parametei^  and  its  density  function  involves  a  matrix 
exponential.  Moment  matching  with  a  nonlinear  programming  approach  is  used  for 
estimating  the  parameters.  This  method  does  not  guarantee  that  the  estimated  cost 
function  will  have  a  unique  and  finite  minimum.  In  this  thesis  we  restric  our  attention  to 
the  case  with  initial  probability  vector  a  =  (1,  0,  0).  Sequential  estimation  using  the 
phase  type  procedure  gives  smaller  average  costs  than  the  nonparametric  procedure  for 
both  small  and  large  sample  sizes  when  the  underlying  long  run  average  cost  function  has 
a  very  distinct  minimum  at  0*.  When  this  cost  function  is  shallow  around  0*,  the 
parametric  procedure  does  not  do  better  than  the  nonparametric  procedure  for  small 
samples.  In  fact,  for  cost  functions  that  are  shallow,  it  is  better  not  to  use  a  maintenance 
policy  i.e.  to  take  $*  =  oo ,  until  the  sample  sizes  are  large  enough  to  give  reasonable 
estimates  of  the  underlying  parameters.  The  reason  for  this  is  that  for  such  cost 
functions,  overestimating  <f>*  does  not  increase  the  long  run  average  cost  much,  at  the 
same  time  the  decrease  in  the  amount  of  censoring  with  over  estimating  0*  greatly 
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improves  estimation  of  the  underlying  distribution  function  F.  Conversely 
underestimating  </>'  for  these  cost  functions  causes  a  drastic  increase  in  long  run  average 
costs  and  increase  censoring  making  estimation  very  difficult. 

Recommendation  for  the  future  research 
The  following  forms  a  list  of  future  research  initiatives: 

-Examine  the  impact  of  the  product-limit  estimator  in  estimating  standardized 
moments  on  the  rest  of  procedure. 

-Do  more  experimentation  on  different  representations  of  phase  type  distributions 
to  get  more  precise  conclusion. 

-Look  for  other  policies  using  phase-sensitive  estimate  and  modified  replacement 
costs  over  time. 

-The  phase-type  sequential  estimations  we  have  been  so  far  assume  the  underlying 
lifetimes  are  iid  and  after  replacement  the  system  is  new.  To  be  more  realistic,  we  can 
modify  the  initial  probability  vector,  increase  transition  rates  between  states,  or  both, 
over  time. 

-Use  phase-type  sequential  estimation  when  the  underlying  life  distribution  F  comes 
from  Gamma,  Weibull  etc.^  that  have  increasing  failure  rate  and  compare  with  the 
nonparametric  procedure.  Since  we  can  use  the  denseness  of  phase  type  distribution 
property  to  approximate  the  set  of  distribution  with  support  on  [0,  a)  i.e.,  the  set  of 
lifetimes,  we  may  get  a  better  method  when  the  distribution  of  lifetimes  is  not  known. 
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APPENDIX  A. 


SHTLE  PH  DISTRIBUTION  ESTIMATION  PARAMETERS 
* 

*  - gams  and  dollar  control  OPTIONS - 

*  (See  Appendices  B  &  C) 

OPTIONS 

WORK  =  100000, 

UMCOL  =  0  ,  UMROW  =  0  ,  SOLPRINT  =  OFF  ,  DECIMALS  =  2 
RESUM  =  100,  ITERIJM  =100000,  OPTCR  =  0.0  ,  STIED  =  3141; 

■*< - DEFINITIONS  AND  DATA- - 

*  This  program  uses  nonlinear  programming  to  estimate  the  3-transient  states  of  phase 

*  type  distribution  parameters.  Matching  the  second  and  third  standardized  moments  are 

*  used.  All  adjustable  data  are  prepared  by  the  "FILE  SCALAR2"  and  uses  the 

*  "SINCLUDE"  statement  bring  to  the  program. 

SINCLUDE  "FILE  SCALAR2" 


—  MODEL 


VARIABLE 

Z  objective  function  value 
A  first  entry  of  the  matrix 
B  second  entry  of  the  matrix 
C  third  entry  of  the  matrix 
D  fourth  entry  of  the  matrix 
E  fifth  entry  of  the  matrix 
F  sixth  entry  of  the  matrix 
G  seventh  entry  of  the  matrix 
H  eighth  entry  of  the  matrix 
I  ninth  entry  of  the  matrix 
ALl  initial  probability  of  being  in  state  1 
AL2  initial  probability  of  being  in  state  2 
AL3  initial  probability  of  being  in  state  3  ; 


EQUATION 

OBJ  Define  objective  function 

EQl  Set  initial  probability  of  being  in  state  1  (nonnegative) 
EQ2  Set  initial  probability  of  being  in  state  2(nonnegative) 
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EQ3  Set  initial  probability  of  being  in  state  3(nonnegative) 

EQ4  Sum  over  all  of  initial  prob.  less  than  or  equal  1.0 

EQS  Set  the  first  diagonal  entry  to  be  nonpositive 

EQ6  Set  the  second  diagonal  entry  to  be  nonpositive 

EQ7  Set  the  third  diagonal  entry  to  be  nonpositive 

EQS  Set  off  diagonal  entries  in  the  first  row  to  be  positive 

EQ9  Set  off  diagonal  entries  in  the  first  row  to  be  positive 

EQIO  Set  sum  off  diagonal  of  the  first  row  less  than  or  equal  first  diagonal 

EQl  1  Set  off  diagonal  entries  in  the  second  row  to  be  positive 

EQ12  Set  off  diagonal  entries  in  the  second  row  to  be  positive 

EQ13  Set  sum  off  diagonal  of  the  second  row  less  than  or  equal  second  diagonal 

EQ14  Set  off  diagonal  entries  in  the  third  row  to  be  positive 

EQ15  Set  off  diagonal  entries  in  the  third  row  to  be  positive 

EQl 6  Set  sum  off  diagonal  of  the  third  row  less  than  or  equal  third  diagonal 

EQ17  Define  the  expected  lifetime  ; 

♦MINIMIZE 

OBJ.. 

Z  =E= 

3*SQR(C  -  C(0)) 

3*P0WER((((2*AL1  *(POWER((-(F*H) +E*I),2) + (C*H-B*I)*(F*G-D*I) 

+ (-(C*E)  +B*F)*(-(E*G) +D*H) 

+ (-(F^H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) + (-(C*E) +B*F)*(B*G-A*H) 

+(-(F*H)+E»I)*(-(C*E)+B»F)+(C*H-B*I)*(C*D-A*F)+(-(C^E)+B*F)*(-(B*D)  + 

A*E)) 

+  2*AL2*((F*G-D*I)*(-(F»H)+E*I)+(-(C*G)+A*I)*(F*G-D*I)+(-(E*G)+D*H)* 
(C*D-A*F) 

+(F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H) 

+(F*G-D*I)*(-(C*E)+B*F)+(-(C*G)+A*I)*(C*D-A*F)+(C*D-A*F)*(-(B*D)+A*E)) 

+  2*AU*((-(E*G)+D*H)*(-(F»H)+E*I)+(B*G-A*H)*(F*G-D*D+(-(B*D)+A*E)* 
(-(E*G)+D*H) 

+ (-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*!) + (-(B*D) + A*E)* 

(B*G-A*H) + (-(E*G) +D*H)*(-(C*E)  +B*F) +(B*G-A*H)*(C*D-A*F)  + 

POWER((-(B*D) + A*E),2)))/ 

POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2) 

-POWER((-(ALl  ♦(-F»H+E*I+C*H-B*I-(C*E) +B*F) 

+AL2*(F*G-D*I-C*G+  A*I+C*D-A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
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(-(C*E*G)  +B*F*G+C*D*H-A*F*H-B*D’*‘1+ A*E*I)),2))**0.5  / 


(-(ALl  *(-(F*H) +E*I+C*H-B*I-(C*E)  +B*F) 

+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 

+AU*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 

(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E’T))  -  MM2), 2) 

***  W2*SQR(GAMMA  -  GAMMA(0)) 

+  P0WER(((((-6*AL1  *((POWER((-(F*H) +E*I),2) + (C*H-B*I)*(F*G-D*I) + 

(-(C*E) +B*F)*(-(E*G) +D*H))*(-(F*H) +E*I) 

+((-(F’^+E*I)*(C*H-B*I)+(C*H-B*I)*(-(C‘^G)+A-T)+(-(C*E)+B*F)*  (B*G-Ani))* 
(F*G-D’"I) 


+ ((-(F*H) +E*I)*(-(C*E)  +B*F) + (C*H-B*D*(C*D-A*F) + (-(C*E)  +B*F)* 

(-(B*D) + A*E))*(-(E*G) +D»H) 

+ (POWER((-(F’^ +E*I),2) +(C-^-B*I)*(F*G-D*I) + (-(C-^ +B*F)-"  (-(E*G) +D*H))* 
(C*H-B*I) 


+ ((-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)’t-(C*G) + A*I) +(-(C*E)  +B*F)*(B*G-A*H))* 
(-(C*G)+A*I) 

+ ((-(F*H) +E*I)*(-(C*E)  +B*F) +(C*H-B*D*(C*D-A*F) +(-(C^E) +B*F)* 

(-(B*D) + A*E))*(B*G-A*H) 


+ (POWER((-(F*H)  +E*I),2) + (C*H-B*I)*(F*G-D*I) + (-(C*E) +B*F)* 
(-(E-^G) +D*H))*(-(C*E)  +B*F) 

+ ((-(F*H) +E’^I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) + (-(C*E) +B*F)* 
(B*G-A*H))*(C*D-A*F) 


+ ((-(F*H) +E*I)*(-(C*E) +B*F) +(C*H-B*I)*(C*D- A*F) + (-(C*E)  +B*F)* 
(-(B’^D) + A*E))*(-(B*D) + A*E)) 

-6*AL2*(((F*G-D*I)*(-(F*H) +E*I) + (-(C*G)+ A*I)*(F*G-D*I) +(-(E*G) +D*H)* 
(C*D-A*F))*(-(F*H) +E*I) 

4-((F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 

(F*G-D*I) 

+ ((F*G-D*I)*(-(C’"E)  +B*F) + (-(C*G) + A*D*(C*D-A'^F) + (C*D-A*F)* 
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(-(B*D) + A*E))*(-(E*G) +D*H) 

+ ((F*G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I) + (-(E*G) +D*H)* 
(C*D-A*F))*(C*H-B*I) 

+ ((F*G-D*I)*(C*H-B*I) +POWER((-(C*G) + A*I)  ,2) + (C*D-A*F)*(B*G-A*H))* 
(-(C*G)+A*I) 

+ ((F*G-D*I)*(-(C*E) +B*F) + (-(C^G) + A*I)*(C*D-A*F) + (C*D-A*F)* 

(-(B*D) + A*E))*(B*G-A*H) 

+ ((F*G-D*I)  *(-(F*H) + E*I) + (-(C*G) + A*I)-"(F*G-D*I) + (-(E*G) + D*H)* 
(C*D-A*F))*(-(C*E) +B*F) 

+ ((F*G-D’"I)*(C-*H-B*I) +POWER((-(C*G)  +  A*I),2) + (C*D-A*F)*(B*G-A*H))* 
(C*D-A*F) 

+ ((F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D- A*F)* 

(-(B*D) + A*E))*(-(B*D) + A*E)) 

-6*AL3*(((-(E*G) +D*H)*(-(F*H) +E*I) +(B*G-A*H)*(F-*G-D*I) + (-(B*D) + A*E)* 
(-(E*G)  +D*H))*(-(F*H) +E*I) 

+ ((-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*I) + (-(B*D)  +  A^E)* 
(B*G-A*H))*(F*G-D*I) 

+((-(E*G) +D*H)*(-(C*E)  +B*F) +(B*G-A*H)*(C*D-A*F)  + 
POWER((-(B*D)+A*E),2))*(-(E*G)+D*H) 

+ ((-(E*G)  +D*H)*(-(F*H) +E*I) + (B*G-A*H)*(F*G-D*I) +(-(B»D) + A*E)* 

(-(E*G) +D*H))*(C*H-B*I) 

+ ((-(E*G) + D*H)*(C*H-B*I) + (B*G- A*H)*(-(C*G) + A*I) + (-(B*D) + A*E)* 
(B*G-A*H))*(-(C*G) + A*I) 

+ ((-(E*G) +D*H)*(-(C*E)  +B*F) + (B*G-A*H)*(C*D-A*F) + 

POWER((-(B*D) + A*E),2))*(B*G-A'^ 

+((-(E*G) +D*H)*(-(F*H) +E*I) +(B*G-A*H)*(F*G-D*I) +(-(B*D) + A*E)* 

(-(E*G) +D*H))*(-(C*E)  +B*F) 

+ ((-(E-^G) +D*H)*(C’^-B*I) + (B*G-A*H)*(-(C*G) + A-T) + (-(B*D) + A*E)* 
(B*G-A*H))*(C*D-A*F) 
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+ ((-(E*G) +D*H)*(-(C*E)  +B*F) + (B*G-A*H)*(C*D-A*F) + 
POWER((-(B*D)+A*E),2))*(-(B*D)+A*E))) 

/POWER((-(C*E*G)  +  B*F*G  +  -  A*F*H  -  B*D*I  +  A*E*I),3) 

+3*((AL1*(-(F*H)+E*H-C*H-B*I-(C*E)+B*F)+AL2*(F*G-D*I-C*G+A*I+C*D- 

A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/(-C*E*G+B*F*G+C*D*H-A*F 

*H-B*D*I+A’"E*I))* 


*♦*  expected  of  lifetime  square 

((2*ALl*(POWER((-(F*H)+E*I),2)+(C*H-B*I)*(F*G-D*I)+(-(C*E)+B*F)* 
(-(E*G) +D*H) + (-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) + 

(-(C*E)  +B*F)*(B*G-A*H) +(-(F*H) +E*I)*(-(C*E)  +B*F) + (C*H-B*I)*(C*D-A*F) 
+ (-(C*E) +B*F)*(-(B*D) + A*E)) 

+2*AL2*((F-"G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I) 

+(-(E*G)  +D*H)*(C*D-A*F) + (F*G-D*I)*(C*H-B*I) +POWER((-(C*G) + A*I),2) 

+ (C*D-A*F)*(B*G-A*H) + (F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) 

+ (C*D-A*F)*(-(B*D) + A*E)) 

+2*AL3*((-(E*G)+D*H)’^(-(F*H)+E*D+(B*G-A*H)*(F’"G-D*I) 

+ (-(B*D) + A*E)*(-(E*G) +D*H) + (-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)* 

(-(C*G) + A*I) +(-(B*D) + A*E)*(B*G-A*H) +(-(E*G) +D*H)*(-(C*E) +B*F) 

+ (B*G-A*H)*(C*D-A*F) +POWER((-(B*D)+ A*E),2)))/ 

POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2)) 


***  2  time  cube  of  the  expected  lifetime  *** 
-2*POWER(((ALl*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F) 

+ AL2*(F*G-D*I-C*G + A*H-  C*D-A*F)+ AL3*(-E*G+D’^+B*G-A*H-B*D+ A*E))/ 
(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I)),3)) 


/((2*AL1  *(POWER((-(F*H)  +E*I),2) + (C*H-B*I)*(F*G-D*I) +(-(C*E)  +B*F)* 
(-(E*G) +D*H) +(-(F*H) +E*I)*(C*H-B*I)+(C’TI-B*I)*(-(C*G) + A*I) + 

(-(C*E)  +B*F)*(B*G-A*H) +(-(F*H) +E*I)*(-(C*E) +B*F) + (C*H-B*I)*(C*D-A*F) 
+ (-(C*E)  +B*F)*(-(B*D) + A*E)) 

+2*AL2*((F*G-D*I)*(-(F*H)+E*I)+(-(C*G)+A*D*(F*G-D*I)+ 

(-(E*G) +D*H)*(C*D-A*F) + (F*G-D*I)*(C*H-B*I) + 

POWER((-(C*G) + A*I),2) + (C*D-A*F)*(B*G-A*H) + (F*G-D*I)* 

(-(C*E) +B*F) + (-(C*G) + A’T)*(C  *D-A*¥) + (C^D- A*F)*(-(B*D) + A*E)) 
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+  2*AL3*((-(E*G) +D*H)*(-(F’Ti) +E*I)+(B*G-A*H)*(F*G-D*I)  + 

(-(B*D)  +  A*E)*(-(E*G)  +  D*H)  +  (-(E*G)  +  D*H)*(C*H-B*I)  +  (B*G-A*H)* 
(-(C*G)  -I-  A*I)  +  (-(B*D)  +  A*E)*(B*G-A»H)  +  (-(E*G) +D*H)*(-(C*E)  +B*F) 
+(B*G-A*H)*(C*D-A*F)+P0WER((-(B*D)+A*E),2)))/ 

(POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2)) 

-P0WER(((AL1  *(-(F*H) +E*I+C*H-B*I-(C*E)  +B*F) 
+AL2*(F*G-D’a-C*G+A*I+C*D-A*F) 

+ AL3*(-E*G +D*H +B*G-A*H-B*D + A*E))/ 

(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I)),2))**1.5)  -  MM3), 2); 
* 

*  SUBJECT  TO 

♦ 

EQl.. 

ALl  =G=  0  ; 

EQ2.. 

AL2  =G=  0  ; 

EQ3.. 

AL3  =G=  0  ; 

EQ4.. 

ALl  +  AL2  +  AL3  =L=  1  ; 

EQ5.. 

A  =L=  0  ; 

EQ6.. 

E  =L=  0  ; 

EQ7.. 

I  =L=  0  ; 

EQ8.. 

B  =G=0; 

EQ9.. 

C  =G=  0 ; 

EQIO.. 

B  +  C  =L=  -A  ; 


EQll.. 

D  =G= 

0; 

EQ12.. 

F  =G= 

0; 

EQ13.. 

D  +  F  =L 

=  -E 

EQ14.. 

G  =G= 

0; 
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EQ15.. 

H  ^G=  0 ; 

EQ16.. 

G  +  H  =L=  -I  ; 

EQ17.. 

-(ALl  *(-(F*H) +E*I + C*H-B*I-(C*E) +B*F) 

+ AU*(F*G-D*I-C*G + A*I + C*D-A*F) 
+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
(-(C*E*G)+B*F*G+C*D*H-A*F*H-B*D*I+A*E*I)  =E=  EX  ; 

- REPORT- - 

MODEL  MATCHM  /ALL/  ; 


*  Parameters  bounding 


A.LO  =  -7.0;  B.LO  =  3.0;  C.LO  =  0.0;  D.LO  =  0.0;  E.LO  =-7.0; 

A.UP  =  -5.0;  B.UP  =  5.0;  C.UP  =  2.0;  D.UP  =  1.0;  E.UP  =-5.0; 

F.LO  =  4.0;  G.LO  =  0.0;  H.LO  =  0.0;  I.LO  =  -8.0; 

F.UP  =  5.0;  G.UP  =  1.0;  H.UP  =  2.0;  I. UP  =  -7.0; 

ALl.LO  =  1.0;  AL2.LO  =  0.0;  AU.LO  =  0.0; 

ALl.UP  =  1.0;  AU.UP  =  0.0;  AL3.UP  =  0.0; 


*  The  initial  solution 


A.L  =  AA  ;  B.L  =  BB  ;  C.L  =  CC;  D.L  =  DD;  E.L  =  EE  ; 
F.L  =  FF  ;  G.L  =  GG  ;  H.L  =  HH  ;  I.L  =  H  ; 

ALl.L  =  ALPl  ;  AL2.L  =  ALP2  ;  AL3.L  =  ALP3  ; 


SOLVE  MATCHM  USING  NLP  MINIMIZING  Z  ; 

DISPLAY  Z.L.A.L  ,  B.L,  C.L,  D.L,  E.L,  F.L,  G.L,  H.L,  I.L, 
ALl.L,  AL2.L,  AL3.L; 

*  Write  output  to  the  CMS  file 

FILE  OUTl  /FILE  OBJ  A  /  ; 

PUT  OUTl  ; 

PUT  Z.L  ; 

FILE  OUT2  /FILE  ALPHA  A  /  ; 

PUTOUT2  ; 

PUT  A.L  ; 

FILE  OUT3  /FILE  BRAVO  A  /  ; 

PUTOUT3  ; 
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PUT  B.L  ; 

FILE  0UT4  /FILE  CHAREE  A  /  ; 
PUT  0UT4  ; 

PUT  C.L  ; 

FILE  OUTS  /FILE  DELTA  A  /  ; 
PUT  OUTS  ; 

PUT  D.L  ; 

FILE  OUT6  /FILE  ECHO  A  /  ; 
PUT  OUT6  ; 

PUT  E.L  ; 

FILE  OUT7  /FILE  FOXTROT  A  / 
PUTOUT7  ; 

PUT  F.L  ; 

FILE  OUTS  /FILE  GOLF  A  /  ; 
PUT  OUTS  ; 

PUT  G.L  ; 

FILE  OUT9  /FILE  HOTEL  A  /  ; 
PUTOUT9  ; 

PUT  H.L  ; 

FILE  OUTIO  /FILE  INDIA  A  /  ; 
PUT  OUTIO  ; 

PUT  LL  ; 

FILE  OUTll  /FILE  ALl  A  /  ; 
PUTOUTll  ; 

PUT  ALl. L  ; 

FILE  OUT12  /FILE  AL2  A  /  ; 
PUT  OUT12  ; 

PUT  AL2.L  ; 

FILE  OUT13  /FILE  AL3  A  /  ; 
PUT  OUT13  ; 

PUT  AU.L  ; 


THE  DnriAL  COMELETE  SYSTEM  LIFETIME  SET  1-5 


N 

SET  1 

SET  2 

SET  3 

SET  4 

SET  5 

1 

0.1226 

0.0740 

0.0855 

0.0528 

0.1143 

2 

0.1455 

0.1413 

0-2157 

0.1139 

0.1999 

3 

0.2662 

0.2261 

0.2209 

0.2412 

0.2002 

4 

0.2925 

0.2629 

0.2839 

0.2892 

0.2422 

5 

0.3205 

0.2831 

0.2852 

0.3118 

0.3164 

6 

0.4214 

0.3340 

0-3116 

0.3485 

0.3398 

7 

0.4394 

0.3911 

0.3403 

0.3973 

0.3534 

8 

0.4685 

0.4617 

0.4517 

0.4442 

0.3893 

9 

0.5804 

0.4761 

0.4691 

0.4686 

0.4410 

10 

0.5871 

0.4909 

0.5254 

0.4777 

0.4766 

11 

0.6228 

0.5290 

0.5536 

0.4885 

0.4964 

12 

0.6534 

0.5818 

0.5993 

0.4970 

0.5074 

13 

0.7220 

0.5962 

0.6575 

0.5114 

0.5240 

14 

0.7590 

0.6233 

0.6623 

0.6073 

0.5825 

15 

0.7618 

0.6419 

0.6855 

0.6433 

0.6033 

16 

0.8206 

0.7407 

0.7156 

0.7576 

0.6255 

17 

0.8245 

0.7682 

0.7690 

0.8057 

0.7391 

18 

0.8317 

0.8312 

0.7866 

0.8740 

0.7723 

19 

1.0056 

0.8693 

0.8852 

1.0004 

0.9128 

20 

1.0259 

0.9727 

0.9473 

1.0424 

1.0832 

21 

1.0853 

1.1084 

1.0478 

1.0634 

1.1936 

22 

1.1659 

1.1272 

1.1155 

1.0659 

1.2003 

23 

1.1852 

1.2037 

1.2936 

1.2238 

1.8976 

24 

1.2040 

1.3319 

1.3906 

1.6404 

2.0769 

25 

1.3069 

2.5362 

1.6715 

3.0033 

2.2631 

E[Z] 

0.705 

0.704 

0.679 

0.735 

0.742 

0.47917 

0.70936 

0.57558 

0.80566 

0.77290 

0.03190 

1.88276 

0.73612 

2.23985 

1.40912 
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THE  INITIAL  COMPLETE  SYSTEM  LIFETIME  SET  6-10 


SET  6 

SET  7 

SET  8 

SET  9 

SET  10 

0.1337 

0.0402 

0.1759 

0.0830 

0.0761 

0.1538 

0.0634 

0.2031 

0.1101 

0.0913 

0.2452 

0.1939 

0.2106 

0.2083 

0.2235 

0.2876 

0.1954 

0.2376 

0.2339 

0.3309 

0.3591 

0.1958 

0.3043 

0.2408 

0.3470 

0.4589 

0.2353 

0.3719 

0.2726 

0.3662 

0.5161 

0.2775 

0.3850 

0.2777 

0.4091 

0.6030 

0.3185 

0.3892 

0.2894 

0.4523 

0.6181 

0.3579 

0.4032 

0.3137 

0.5581 

0.6761 

0.3922 

0.4199 

0.3227 

0.5800 

0.6851 

0.3985 

0.6029 

0.3390 

0.5805 

0.7738 

0.5206 

0.6406 

0.3873 

0.5993 

0.8541 

0.5233 

0.6506 

0.3878 

0.6354 

0.8654 

0.6261 

0.7156 

0.4055 

0.6571 

0.9539 

0.6336 

0.7584 

0.4076 

0.7529 

1.0193 

0.6776 

0.7677 

0.4333 

0.7874 

1.0352 

0.8705 

0.8754 

0.6178 

0.7991 

1.0677 

1.1124 

0.8836 

0.6245 

0.8572 

1.1196 

1.1228 

0.9575 

0.6455 

0.8920 

1.1995 

1.1742 

0.9805 

0.9824 

0.9024 

1.2301 

1.1872 

1.3367 

1.0054 

0.9699 

1.3043 

1.1887 

1.4040 

1.1143 

0.9734 

2.0273 

1.4849 

1.5844 

1.2279 

1.0452 

3.2891 

1.7218 

1.6250 

1.2420 

1.1461 

3.5360 

2.5301 

1.8392 

2.1319 

1.1594 

1.000 

0.722 

0.749 

0.572 

0.648 

0.82457 

0.81263 

0.63261 

0.80948 

0.46844 

1.89734 

1.25687 

0.79102 

1.69217 

0.14707 

8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


E[Z] 


2“* 


3 
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INITIAL  COMPLETE  SYSTEM  LIFETIME  SET  11-15 


N 


8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


E[Z] 


2nd 


SET  11 

SET  12 

SET  13 

SET  14 

SET  15 

0.0268 

0.0433 

0.0470 

0.1730 

0.1101 

0.0745 

0.1285 

0.0922 

0.2556 

0.1386 

0.1235 

0.2052 

0.1295 

0.3124 

0.1709 

0.1800 

0.2435 

0.2064 

0.3125 

0.1998 

0.2526 

0.2924 

0.2924 

0.3447 

0.2075 

0.2704 

0.3335 

0.2957 

0.3894 

0.2588 

0.2972 

0.3645 

0.3110 

0.4044 

0.2669 

0.3697 

0.4936 

0.3209 

0.5555 

0.2737 

0.3811 

0.4984 

0.4173 

0.5796 

0.2905 

0.3990 

0.5064 

0.4577 

0.5803 

0.3506 

0.4070 

0.5778 

0.5168 

0.5928 

0.3604 

0.4855 

0.6699 

0.5179 

0.6056 

0.4117 

0.6402 

0.7442 

0.5258 

0.6525 

0.5809 

0.7087 

0.7502 

0.5456 

0.7299 

0.5981 

0.7408 

0.8656 

0.5547 

0.7355 

0.6558 

0.8496 

0.8831 

0.5559 

0.8122 

0.7306 

0.9290 

1.0607 

0.5633 

0.9379 

0.7388 

0.9873 

1.1003 

0.5692 

0.9433 

0.7400 

0.9983 

1.1834 

0.6501 

1.0197 

0.7786 

1.0318 

1.2235 

0.7978 

1.2675 

0.7936 

1.0937 

1.3370 

0.8750 

1.3842 

1.1532 

1.1137 

1.5754 

1.0843 

1.6525 

1.2444 

1.2706 

1.8085 

1.1683 

1.8276 

1.2781 

1.4875 

1.9023 

1.1881 

1.8545 

1.3428 

1.4958 

1.9894 

1.3242 

2.8764 

1.5937 

0.665 

0.831 

0.560 

0.872 

0.611 

0.64986 

0.66778 

0.60540 

0.71398 

0.68553 

0.32858 

0.60712 

0.69085 

1.51250 

0.77568 
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THE  DJITIAL  COMPLETE  SYSTEM 


N 


8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


E[X] 


2“* 


3 


SET  16 

SET  17 

SET  18 

SET  19 

SET  20 

0.1677 

0.0565 

0.0398 

0.1798 

0.0202 

0.1819 

0.1811 

0.1402 

0.2037 

0.1920 

0.2139 

0.1941 

0.1506 

0.2197 

0.2205 

0.3436 

0.2810 

0.1792 

0.3159 

0.2209 

0.3641 

0.3315 

0.1922 

0.3464 

0.3406 

0.4393 

0.4147 

0.3787 

0.3629 

0.3645 

0.4800 

0.4249 

0.4660 

0.3744 

0.3941 

0.5328 

0.4629 

0.4944 

0.3793 

0.4130 

0.5643 

0.4939 

0.5952 

0.3870 

0.5525 

0.7220 

0.5382 

0.6198 

0.3875 

0.5729 

0.7393 

0.5776 

0.6269 

0.4011 

0.5787 

0.7666 

0.6003 

0.6276 

0.4906 

0.5835 

0.7766 

0.6055 

0.6608 

0.7436 

0.6138 

0.7964 

0.6230 

0.7183 

0.7690 

0.6389 

0.7981 

0.6477 

0.7970 

0.7824 

0.6609 

0.8228 

0.7631 

0.7997 

0.8226 

0.6624 

0.8390 

0.7980 

0.8065 

0.8492 

0.7091 

0.9422 

0.7988 

0.8065 

0.9542 

0.7178 

1.0115 

0.8518 

0.9682 

1.2333 

0.7447 

1.0381 

0.9091 

0.9939 

1.3601 

1.2027 

1.1195 

1.0447 

1.0028 

1.4489 

1.3247 

1.2716 

1.1973 

1.0875 

1.5989 

1.4412 

1.3426 

1.2331 

1.2722 

2.2702 

1.6155 

1.6082 

1.2339 

1.4913 

2.2971 

1.7162 

1.7724 

1.3386 

1.8663 

3.1519 

2.1912 

0.786 

0.664 

0.711 

0.893 

0.748 

0.52328 

0.51958 

0.60066 

0.83695 

0.70314 

0.55452 

0.32242 

0.66254 

1.46868 

1.13073 
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APPENDK  B. 

PROGRAM  COST 
C 

C  THIS  PROGRAM  PROVIDES  THE  OPTIMAL  AGE  REPLACinviLNT  OF  THE 
C  PHASE  TYPE  DISTRIBUTION.  IT  WORKS  ONLY  SPECIFIC  CASE  OF  3 
C  PHASESOFTRANSIENTSTATES.  THE  PARAMETERS  ARE  AS  FOLLOWS: 
C  T(I,J)  =  ENTRIES  OF  THE  TRANSITION  MATRIX  OF  PH  DISTN. 

C  AL(I)  =  INITIAL  PROBABILITY  OF  BEING  IN  STATE  I 

C  Cl  =  PLANNED  MAINTENANCE  COST 

C  C2  =  UNPLANNED  MAINTENANCE  COST 

C  L  =  MAXIMUM  LIFETIME  THAT  WANT  TO  CALCULATE 

C  DET  =  DETERMINANT  OF  TRANSITION  MATRIX 

C  SMALL  =  THE  OPTIMAL  AGE  REPLACEMENT 

C  CS  =  AGE  REPLACEMENT  COSTS 
* 

*  INITIAL  INPUT 

•K 

INTEGER  I,J,K,N 

REAL  Cl,  C2,  TI(3,3),  T(3,3),  DET,  AL(3),  L,  OBJ 
REAL  TS(3,3),  SMALL,  PI,  A,  IDEN(3,3),  S,  NUMER 
REAL  DENOM,  P(3,3),  Q(3,3),  CS(400),  EXPTS(3,3) 

DATA  C1,C2,L,TS/100.0,500.0,4.0,9*0.0  / 

N  =  0 
S  =  0.01 

OPEN(UNIT  =11, FILE  =  ’ALPHA’) 

OPEN(UNIT  =12,FILE  =  ’BRAVO’) 

OPEN(UNIT  =13,FILE  =  ’CHAREE’) 

OPEN(UNIT  =14,FILE  =  ’DELTA’) 

OPEN(UNIT  =15,FILE  =  ’ECHO’) 

OPEN(UNIT  =16,FILE  =  ’FOXTROT’) 

OPEN(UNrr  =17,FILE  =  ’GOLF’) 

OPEN(UNIT  =18, FILE  =  ’HOTEL’) 

OPEN(UNIT  =19,FILE  =  ’INDIA’) 

OPEN(UNIT  =20,FILE  =  ’ALl’) 

OPEN(UNIT  =21, FILE  =  ’AL2’) 

OPFJ4(UNIT  =22,FILE  =  ’AL3’) 

OPEN(UNIT  =23,FILE  =  ’OBJ’) 

READ(11,*)T(1,1) 

READ(12,*)T(1,2) 

READ(13,*)T(1,3) 


47 


READ(14,*)T(2,1) 

READ(15,*)T(2,2) 

READ(16,*)T(2,3) 

READ(17,*)T(3,1) 

READ(18,*)T(3,2) 

READ(19,*)T(3,3) 

READ(20,*)AL(1) 

READ(21,*)AL(2) 

READ(22,*)AL(3) 

READ(23,*)OBJ 

IF(OBJ  .GT.  0.05)  GOTO  555 

*  CALCULAnON  COST  FOR  EACH  REPLACEMENT  TIME 

•K 

DET  =  -(T(1,3)’^(2,2)*T(3,1)) 

/  +  T(l,2)*T(2,3)n'(3,l) 

/  +  T(1,3)*T(2,1)’T(3,2) 

/  -T(l,l)*T(2,3)n'(3,2) 

/  -T(1,2)*T(2,1)*T(3,3) 

/  +  T(1,1)-"T(2,2)*T(3,3) 

PRINT*, ’DET  =  ’,DET 

TI(1,1)  =  (-(T(2,3)*T(3,2))  +  T(2,2)*T(3,3))/DEr 
TI(1,2)  =  (T(1,3)*T(3,2)-T(1,2)*T(3,3))/DET 
TI(1,3)  =  (-(T(1,3)*T(2,2))  +  T(1,2)*T(2,3))/DET 
TI(2,1)  =  (T(2,3)*T(3,1)  -  T(2,1)*T(3,3))/DET 
TI(2,2)  =  (-(T(1,3)*T(3,1))  +  T(1,1)*T(3,3))/DET 
11(2,3)  =  (T(1,3)*T(2,1)  -  T(1,1)*T(2,3))/DET 
11(3,1)  =  (-(T(2,2)*T(3,1))  +  T(2,1)*T(3,2))/DET 
11(3,2)  =  (T(1,2)*T(3,1)-T(1,1)*T(3,2))/DET 
11(3,3)  =  (-(T(1,2)*T(2,1))  +  T(1,1)*T(2,2))/DET 
DO  101  =  1,3 
DO  5  J  =  1,3 

IF(I.EQ.J)THEN 
IDENa,J)  =  1.0 
ELSE 

IDENa,J)  =  0.0 
ENDIF 

5  CONTINUE 
10  CONTINUE 

♦ 

100  CONTINUE 
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N  =  N+1 
DO  20  I  =  1,3 
DO  15  J  =  1,3 
TSa.J)  =  S*Ta,J) 

15  CONTINUE 

20  CONTINUE 

* 

*  FIND  THE  EXPONENTIAL  OF  MATRIX  TS 

* 

CALL  EXPON(TS,EXPTS,PI) 

A  =  0.0 
DO  21  I  =  1,3 
D0  22  J  =  1,3 

A  =  A  +  AL(I)*EXPTS(I,J) 

22  CONTINUE 

21  CONTINUE 

NUMER  =  C2  -  A*(C2-C1) 

DO  30  I  =  1,3 
DO  25  J  =  1,3 

P(I,J)  =  EXPTS(I,J)  -  IDEN(I,J) 

25  CONTINUE 
30  CONTINUE 
DO  45  I  =  1,3 
DO40  J  =  1,3 
Qa,J)  =  0.0 
D0  35K  =  1,3 

Qa,j)  =  Qa,J)  +Tia,K)*p(K,j) 
35  CONTINUE 
40  CONTINUE 
45  CONTINUE 
DENOM  =  0.0 
DO  55  I  =  1,3 
DO50  J  =  1,3 

DENOM  =  DENOM  +  AL(I)*Q(I,J) 
50  CONTINUE 
55  CONTINUE 

CS(N)  =  NUMER/DENOM 
S  =  S+0.01 
IF(S.LE.L)GOTO  100 

*  FIND  THE  TIME  THAT  HAS  MINIMUM  COST 
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K  =  1 

200  IF(K.EQ.N)THEN 
SMALL  =  K*0.01 
ELSEIF(CS(K+ 1).GT.CS(K))THEN 
SMALL  =  K*0.01 
ELSEIF(CS(K+  1).LE.CS(K))TEIEN 
SMALL  =  (K+l)*0.01 
K  =  K  +  1 
GOTO  200 
ENDDF 

♦ 

*  OUTPUT  OPTIMAL  REPLACEMENT  AGE 

* 

OPEN(UNIT  =  3, FILE  =  ’OPTAGE’) 

WRITE(3,333)SMALL  ,  PI 
333  FORMAT(lX,F6.3,4X,E11.4) 

OPEN(UNIT  =  4, FILE  =  ’OPTCOST’) 

WRITE(4,*)CS(K) 

555  STOP 
END 

SUBROUTINE  EXPON(A,EXPA,PI) 

THIS  SUBROUTINE  USES  TO  FIND  THE  EXPONENTIAL  OF  A  3X3  MATRIX. 
IT  WORK  WITH  COUPLE  OF  IMSL  SUBROUTINES 
(’EVCRG’ ,  ’UNCG’ ,  ’EPIRG’). 

INTEGER  LDA,LDEVEC,N,NOUT,I,J,K 

PARAMETER(N=3,LDA=N,LDEVEC=N,LDUINV=N,LD=N,LDU=N) 
REAL  A(LDA,N),PI,EXPA(3,3) 

COMPLEX  EVAL(N),  EVEC(LDEVEC,N),  ELD(3,3),  UINV(3,3), 
COMPLEX  C(3,3),  D(3,3) 

OPEN(UNIT=  1  ,FILE=  ’EXPONENT’) 

*  CALCULATE  THE  EIGENVALUES  AND  EIGENVECTERS 

* 

CALL  EVCRG(N,A,LDA,EVAL,EVEC,LDEVEC) 

DO  15  I  =  1,3 
DO  10  J  =  1,3 
ELD(I,J)  =  0.0 
10  CONTINUE 
15  CONTINUE 
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ELD(1,1)  =  CEXP(EVAL(1)) 

ELD(2,2)  =  CEXP(EVAL(2)) 

ELD(3,3)  =  CEXP(EVAL(3)) 

* 

*  CALCULATE  THE  INVERSE  EIGENVECTERS’  MATRIX 

iK 

CALL  LINCG(N,EVEC,LDU,UINV,LDUINV) 

DO  30  I  =  1,3 
D0  25  J  =  1,3 
Ca,J)  =  0.0 
DO  20  K  =  1,3 

C(I,J)  =  C(I,J)  +  ELD(I,K)*UINV(K,J) 

20  CONTINUE 

25  CONTINUE 
30  CONTINUE 
DO  45  I  =  1,3 
DO40  J  =  1,3 
D(I,J)  =  0.0 
DO  35  K  =  1,3 

D(I,J)  =  D(I,J)  +  EVEC(I,K)»C(K,J) 

35  CONTINUE 

40  CONTINUE 
45  CONTINUE 
DO  55  I  =  1,3 
DO50  J  =  1,3 

EXPAa.J)  =  REAL(D(I,J)) 

50  CONTINUE 
55  CONTINUE 
C 

C  CALCULATE  THE  PERFORMANCE  INDEX  OF  EIGENVALUES  AND 
C  EIGENVECTERS  IF  LESS  THAN  1  ’EXCELLENT’, IF  IT  IS  BETWEEN  1 
C  AND  100  ’GOOD’,  IF  ITIS  GREATER  THAN  100  ’THE  SOLUTION  IS 
C  SUSPECTED’ 

C 

PI  =  EPIRG(N,N,A,LDA,EVAL,EVEC,LDEVEC) 

RETURN 

END 
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PROGRAM  SIMUF 


C 

C  TfflS  PROGRAM  PROVIDES  A  RANDOM  NUMBER  OF  A  PHASE  TYPE 
C  DISTRIBUTION  THAT  HAVE  SPEOTIC  REPRESENTATION  WITH  AN 
C  INITIAL  PROBABHITIES  VECTOR(ALPHA)  AND  A  MATRIX  OF 
C  TRANSITION  RATES  BETWEEN  TRANSIENT  STATES.  THIS  PROGRAMS 
C  WORKS  WITH  SUBROUTINE  ’RANNUM’  AND  MORE  SPECIFIC  CASE 
C  ONLY  4  STATES  DISTRIBUTION.  THE  VARIABLES  ARE  AS  FOLLOWS: 

C  LD(I, J)  =  TRANSITION  RATE  FROM  STATE  I  TO  STATE  J 

C  P(I,J)  =  TRANSITION  PROBABUJTY  FROM  STATE  I  TO  STATE  J 
C  AL(I)  =  INITIAL  PROBABILITY  OF  BEING  IN  STATE  I 

C  EUF  =  THE  SUPPOSED  UPPER  BOUND  EXPECTED  LIFETIME 

* 

*  INITIALIZATION 

He 

INTEGER  S,SEED,N,N1 
REAL  U,EXP,LIF,ELIF 

REALLD12,LD13,LD14,LD21,LD23,LD24,LD31,LD32,LD34 
REAL  P12,P13,P21,P23,P31,P32,AL1,AL2,AL3,D 
DATA  LD12,LD13,LD14,LD21  ,LD23,LD24,LD31  ,LD32,LD34 
+  /5.0,5.0,5.0,6.0,6.0,6.0,7.0,7.0,7.0/ 

DATA  P12,P13,P21,P23,P31,P32,AL1,AL2,AU, EUF 
+  70.8,0.2,0. 167,0.833,0. 143,0.286, 1 .0,0.0,0.0, 10.0/ 

CALL  EXCMS(’FILEDEF  10  DISK  FILE  LFTEST  (DISP  MOD’) 
OPEN(UNIT  =  1,FILE  = ’LIFETIME’) 

OPEN(UNIT  =  2,FILE  =’SEED’) 

READ(2,*)SEED 

CLOSE(2) 

UF  =  0.0 


GENERATE  UNIFORM  0  AND  1  FOR  INITIAL  PROBABILITY 

* 

CALL  RANNUM(1,SEED,0.0,1.0,0,U) 

IF  (U.LE.AL1)  THEN 
S  =  1 

ELSEIF(U.LE.  ALl  +  AL2)THEN 
S  =  2 
ELSE 
S  =  3 
ENDIF 

* 

*  SIMULATE  LIFETIME  WHEN  IS  IN  STATE  1 
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5  CONTINUE 

IF  ((S.EQ.l)  .AND.  (UF.LT.EUF))  THEN 

CALL  RANNUM(1, SEED, 0.0,1. 0,0, U) 
IF(U.GT.P12+P13)THEN 
S  =  4 

CALL  RANNUM(3,SEED,LD14,0,0,EXP) 
UF  =  LIF  +  EXP 
ELSEIF(U.GT.P12)THEN 
S  =  3 

CALL  RANNUM(3,SEED,LD13,0,0,EXP) 
UF  =  UF  +  EXP 
ELSE 
S  =  2 

CALL  RANNUM(3,SEED,LD12,0,0,EXP) 
UF  =  UF  +  EXP 
ENDIF 
ENDIF 

* 

*  SIMULATE  UFETIME  WHEN  IS  IN  STATE  2 

* 

IF((S.EQ.2)  .AND.  (UF.LT.ELIF))THEN 
CALL  RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P21  +P23)THEN 
S  =  4 

CALL  RANNUM(3,SEED,LD24,0,0,EXP) 
UF  =  UF  +  EXP 
ELSEIF(U.GT.P21)THEN 
S  =  3 

CALL  RANNUM(3,SEED,LD23,0,0,EXP) 
UF  =  UF  +  EXP 
ELSE 
S  =  1 

CALL  RANNUM(3,SEED,LD21,0,0,EXP) 
UF  =  UF  +  EXP 
ENDIF 
ENDIF 

>*< 

*  SIMULATION  LIFETIME  WHEN  IS  IN  STATE  3 

♦ 

IF((S.EQ.3)  .AND.  (LIF.LT.ELIF))THEN 
CALL  RANNUM(1, SEED, 0.0, 1.0,0,U) 
IF(U.GT.P31  +P32)THEN 
S  =  4 
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CALL  RANNUM(3, SEED, LD34, 0,0, EXP) 

UF  =  UF  +  EXP 
ELSEIF(U.GT.P31)THEN 
S  =  2 

CALL  RANNUM(3,SEED,LD32,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSE 
S  =  1 

CALL  RANNUM(3,SEED,LD31,0,0,EXP) 

TTF  =  UF  +  EXP 
ENDIF 
ENDIF 

IF((S.LT.4)  .AND.  (LIF.LT.ELIF))GOTO  5 

* 

*  OUTPUT  TO  ’FILE  LIFETIME  A’ 

WRrTE(10,*)LIF 

WRITE(1,*)LIF 

OPEN(UNIT  =  2,FILE  =’SEED’) 

WRITE(2,*)SEED 

STOP 

END 

nr «  **«««««*«****  *4>  **«««««  «*«««**«  lliltl « W  *  4t  *««*«*«  Dull «  Hr  *«  *****  **  >l> 

SUBROUTINE  RANNUM(DISTN,SEED,RPARM1,RPARM2,IPARM,X) 

C 

C  THIS  SUBROUTINE  PROVIDES  A  RANDOM  NUMBER  OF  7  DISTRIBUTION. 
C  IT  INTERFACES  WITH  THE  LLRANDOMB  ROUTINES  PROVIDED  IN  THE 
C  NONIMSL  LIBRARY.  THE  PARAMETER  AND  CALUNG  PROCEDURE  ARE 
C  AS  FOLLOWS: 

C  DIST  =  DISTRIBUTION  TYPE  YOU  WANT  TO  SELECT  AN  INTEGER 
C  BETWEEN  1  AND  7. 

C  SEED  =  THE  RANDOM  NUMBER  SEED  YOU  WISH  TO  USE. 

C  RPARM1,RPARM2,  AND  IPARM  ARE  REAL  AND  INTEGER  PARAMETERS. 
C  PASSED  TO  THE  ROUTINE  WITH  MEANINGS  WHICH  VARY  WITH  THE 
C  TYPE  OF  DISTRIBUTION  YOU  DESIRE. 

C  X  =  THE  RETURNED  RANDOM  NUMBER,  IT  IS  ALWAYS  REAL. 

C  DISTRIBUTION  NUMBERS  AND  THE  ASSOCIATED  FARM  DEFINITIONS 
C  1  =  UNIFORM  ON  THE  INTERVAL  RPARMl  TO  RPARM2 
C  2  =  NORMAL  WITH  MEAN  RPARMl  AND  VARIANCE  RPARM2 
C  3  =  EXPONENTIAL  WITH  RATE  RPARMl 
C  4  =  COUCHY  WITH  A  =  RPARMl  AND  B  =  RPARM2 
C  5  =  GAMMA  WITH  SHAPE  RPARM2  AND  RATE  RPARMl 
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C  6  =  POISSON  WITH  RATE  RPARMl 
C  7  =  GEOMETRIC  WITH  P  =  RPARMl 
C 

REAL  RPARMl, RPARM2,X 
INTEGER  DISTN,SEED,IPARM,N 
REAL  TEMP,VARIAT(1) 

IF(DISTN.LE.O  .OR.  DISTN.GE.8)THEN 
WRITE(10,'")’DISTN  DOES  NOT  EXIT’ 

STOP 

ENDIF 

* 

GOTO  (10,20,30,40,50,60,70),DISTN 

* 

*  GENERATE  A  UNIFORM  BETWEEN  RPARMl  AND  RPARM2 
10  CONTINUE 

IF(RPARM1-RPARM2.EQ.0)THEN 
WRITE(10,*)’ILLEGAL  EQUAL  RPARMS  IN  RANNUM’ 

STOP 

ENDIF 

IF(RPARM1  .GT.RPARM2)THEN 
TEMP  =  RPARMl 
RPARMl  =  RPARM2 
RPARM2  =  TEMP 
ENDIF 

CALL  LRND(SEED,  VARIAT,  1,1,0) 

VARIAT(l)  =  RPARMl  +  (RPARM2-RPARM1)*VARIAT(1) 

GOTO  99 

*  GENERATE  A  NORMAL  WITH  MEAN  RPARMl  AND  STD.  DEV.  RPARM2 

* 

20  CALL  LNORM(SEED,  VARIAT,  1,1,0) 

WRITE(10,*)’NORMAL(0,1)  ’,VARIAT(1) 

VARIAT(l)  =  (VARIAT(1)*RPARM2)  +  RPARMl 
GOTO  99 

* 

*  GENERATE  AN  EXPONENTIAL  WITH  RATE  RPARMl  (1/MEAN) 

30  CONTINUE 

IF(RPARM1  .EQ.0.0)THEN 
WRITE(10,*)’ILLEGAL  ZERO  RATE  IN  RANNUM’ 

STOP 

ENDIF 
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CALL  LEXPN(SEED,VARIAT,  1,1,0) 

VARIAT(l)  =  VARIAT(1)/RPARM1 
GOTO  99 

GENERATE  A  COUCHY  WITH  A  =  RPARMl  AND  B  =  RPARM2 

40  CONTINUE 

IF(RPARM2.LE.0.0)THEN 

WRITE(10,*)’ILLEGAL  COUCHY  SPREAD  IN  RANNUM,B=  ’,RPARM2 
STOP 
ENDIF 

CALL  LCCHY(SEED,VARIAT,  1,1,0) 

VARIAT(l)  =  (VARIAT(1)*RPAEM2)  +  RPARMl 
GOTO  99 

GENERATE  A  GAMMA  WITH  SHAPE  RPARM2  AND  RATE  RPARMl 

50  CONTINUE 

IF(RPARM1  .LE.0.0)THEN 

WRITE(10,*) ’ILLEGAL  NONPOSmVE  GAMMA  RATE  IN  RANNUM’ 
STOP 
ENDIF 

IF(RPARM2.LE.0.0)THEN 

WRITE(10,*)’ILLEGAL  SHAPE  PARAMEDER  IN  RANNUM’ 

STOP 

ENDIF 

CALL  LGAMA(SEED,VARIAT,1,1,0,RPARM2) 

VARIAT(l)  =  VARIAT(1)*(1.0/RPARM1) 

GOTO  99 

GENERATE  A  POISSON  WITH  RATE  RPARMl 

60  CONTINUE 

IF(RPARM1  .LE.0.0)THEN 

WRITE(10,*)’ILLEGAL  POISSON  RATE  IN  RANNUM’ 

STOP 

ENDIF 

CALL  LPOIS(SEED,VARIAT,1,1,0,RPARM1) 

GOTO  99 

GENERATE  A  GEOMETRIC  WITH  P  =  RPARMl 


70  CONTINUE 


IF(RPARM1  .LE.O.O)THEN 

WRITE(10,*)’ILLEGAL  GEOM.  PROB.  IN  RANNUM’ 
STOP 
ENDIF 

CALL  LGEOM(SEED,VARIAT,1,1,0,RPARM1) 

GOTO  99 

99  CONTINUE 
X  =  VARIAT(l) 

RETURN 

END 
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PROGRAM  MOMENT 


THIS  PROGRAM  USES  TO  COMPUTE  THE  SECOND  AND  THIRD 
STANDARDIZED  MOMENTS  THAT  SUPPORTS  MOMENT  MATCHING  IN 
GAMS  PROGRAM. 

•  INITIAUZE  DATA 

* 

INTEGER  UMIT,N,Nl,I,J,COUNT 

REAL  UF,  NEW(2),  REPAGE,  A,  B,  C,  EX,  EXSQ,  EXC,  EX3 
REAL  STD,  SECOND,  THIRD,  Cl,  C2,  AAC,  OBJ 
REAL  RAWDAT(1000,2),  FB(IOOO),  P(IOOO),  TOT,  TCOST 
DATA  Cl  ,C2,TOT,TCOST, AAC/100.0,500.0,0.0,0.0,0.0/ 

OPEN(UNIT  =  1,FILE  =’ LIFETIME’,  STATUS  =’OLD’) 

OPEN(UNIT  =  2,FILE  = ’OPTAGE’ , STATUS  =’OLD’) 

OPEN(UNIT  =  3, FILE  =’RAWDATA’,  STATUS  =’OLD’) 

OPEN(UNIT  =  4,FILE  = ’MOMENT’) 

OPEN(UNIT  =  5,FILE  =’ COUNT’,  STATUS  =’OLD’) 

READ(1,’*‘)LIF 
READ(2,*)REPAGE 
READ(5,*)COUNT 
PRINT*, ’COUNT= ’, COUNT 
N1  =  25  +  COUNT 
DO  5  I  =  1,1000 
P(I)  =  0.0 
FB(I)  =  0.0 
5  CONTINUE 
DO  44  I  =  1,N1-1 

READ(3, 1 1)RAWDAT(I,  1),RAWDAT(I,2) 

11  F0RMAT(1X,F7.4,4X,F3.1) 

44  CONTINUE 
CLOSE(3) 

OPEN(UNIT=3,FILE=  ’RAWDATA’) 

« 

*  CALCULATION 

IF(LIF.LE.REPAGE)THEN 
NEW(l)  =  LIF 
NEW(2)  =  1.0 
ELSE 

NEW(l)  =  REPAGE 
NEW(2)  =  0.0 
ENDIF 
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CALL  INSERT(N1,NEW,RAWDAT) 

DO  55  I  =  1,N1 

WRITE(3,33)  RAWDATa,l),RAWDATa,2) 

33  F0RMAT(1X,F7.4,4X,F3.1) 

55  CONTINUE 

IF(RAWDAT(1,2)  .NE.  0.0)THEN 
FB(1)  =  (N1-1.0)/N1 
ELSE 

FB(1)  =  1.0 
ENDIF 

DO  10  I  =  2,N1 
IF(RAWDATa,2).EQ.0.0)THEN 
FB(I)  =  FB(I-1) 

ELSEIF(I.LT.N1)THEN 

FB(I)  =  ((Nl-I)/(Nl  +  1.0-I))*FBa-l) 

ELSE 

FB(N1)  =  0.0 
ENDIF 

10  CONTINUE 

P(l)  =  l.O-FB(l) 

DO  22  I  =  2,N1 
P(I)  =  FBa-l)-FB(I) 

22  CONTINUE 
EX  =  0.0 
EXSQ  =  0.0 
EXC  =  0.0 
EX3  =  0.0 
DO  15  J  =  1,N1 
TOT  =  TOT  +  RAWDAT(J,1) 

TCOST  =  TCOST  +  C2*RAWDAT(J,2)  +  C1^(1-RAWDAT(J,2)) 
A  =  RAWDAT(J,1)*P(J) 

B  =  (RAWDAT(J,1)**2)*P(J) 

C  =  (RAWDAT(J,1)**3)*P(J) 

EX  =  EX  +  A 
EXSQ  =  EXSQ  +  B 
EXC  =  EXC  +  C 
15  CONTINUE 

AAC  =  TCOST/TOT 

STD  =  (EXSQ  -  EX**2)**0.5 

EX3  =  EXC  +  (2*EX**3)  -  (3*EX*EXSQ) 

SECOND  =  STD/EX 
THIRD  =  EX3/(STD**3) 
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*  OUTPUT 

* 

WRITE(4, 100)SECOND, THIRD,  AAC, EX 
100  FORMAT(2X,Fl  1 .8,4X,F1 1 . 8,4X,F8.3,4X,F6.3) 

*  CLOSE(5) 

*  OPEN(UNIT  =  5, FILE  = ’COUNT’, STATUS  =’OLD’) 

*  WRrrE(5,*)COUNT+l 
STOP 

END 

iifitt*>t‘*******************************4‘*#*****>t‘***********)t‘************ 

SUBROUTINE  INSERT(LIMIT,NEW,X) 

C  THIS  SUBROUTINE  USES  TO  INSERT  A  PAIR  ’NEW’  DATA  INTO  THE 
C  DATA  FILE  ’X’  IN  ASCENDING  ORDER. 

INTEGER  LIMIT,J,K 
REAL  NEW(2),X(1000,2) 

LOGICAL  DONE 
DONE  =  .FALSE. 

J  =  1 

5  IF((J.LT.LIMIT).AND.(,NOT.DONE))THEN 
IF(X(J,  1).LT.NEW(1))THEN 
J  =J+1 
ELSE 

DONE  =  .TRUE. 

ENDIF 
GOTO  5 
ENDIF 

IF(J.EQ.LIMIT)THEN 
X(J,1)  =  NEW(l) 

X(J,2)  =  NEW(2) 

ELSE 

IF(J.LT.LIMIT)THEN 
DO  lOK  =  LIMIT,J+1,-1 
X(K,1)  =  X(K-1,1) 

X(K,2)  =  X(K-1,2) 

10  CONTINUE 

X(J,1)  =  NEW(l) 

X(J,2)  =  NEW(2) 

ENDIF 

ENDIF 

RETURN 

END 
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PROGRAM  DATA 


C 

C  TfflS  PROGRAM  PROVIDES  THE  RESULT  SEQUENTIAL  ESTIMATED 
C  ENTRIES  TRANSITION  RATES  Ta,J)  AND  INITIAL  STATE  PROB. 

C  ALPAH(I),AND  ALSO  ADJUST,  PREPARE  ALL  THE  INITIAL  DATA  THAT 
C  WILL  BE  USED  FOR  THE  NEXT  SEQUENTIAL  ESTIMATION. 

C  THE  VARIABLE  ARE  AS  FOLLOWS: 

C  A,B,C 

C  D,E,F  =  THE  ESTIMATED  TRANSITION  MATRIX 
C  G,H,I 

C  AL1,AL2,AL3  =  THE  ESTIMATED  INITIAL  STATES  PROBABIUnES 
C  X  =  RANDOM  VARIABLE  OF  LIFETIME 
C  OPTAGE  =  OPTIMAL  AGE  REPLACEMENT 
C  SECOND  =  THE  SECOND  STANDARDIZED  MOMENT 
C  THIRD  =  THE  THIRD  STANDARDIZED  MOMENT 
C  OBJ  =  OBJECTIVE  FUNCTION  VALUE  OF  GAMS  PROGRAM 

■K 

*  INPUT 

* 

INTEGER  COUNT 

REAL  A,B,C,D,E,F,G,H,I,AL1,AL2,AU,X 

REAL  EX, OPTCS, OPTAGE, Z,SECOND,THIRD,DI, OBJ, AAC, PI 

OPEN(UNIT  =  10, FILE  =’OBJ’,STATUS  =’OLD’) 

OPENCUNTT  =  11, FILE  =’ ALPHA’, STATUS  =’OLD’) 

OPEN(UNIT  =  12,FILE  = ’BRA  VO’,  STATUS  =’OLD’) 

OPEN(UNIT  =  13, FILE  =’CHAREE’,  STATUS  =’OLD’) 

OPENCUNTT  =  14,FILE  = ’DELTA’, STATUS  =’OLD’) 

OPENCUNTT  =  15,FILE  = ’ECHO’, STATUS  =’OLD’) 

OPENCUNTT  =  16,FILE  =’ FOXTROT’,  STATUS  =’OLD’) 

OPENCUNTT  =  17,FILE  = ’GOLF’, STATUS  =’OLD’) 

OPENCUNTT  =  18,FILE  = ’HOTEL’, STATUS  =’OLD’) 

OPENCUNTT  =  19, FILE  = ’INDIA’, STATUS  =’OLD’) 

OPENCUNTT  =  20,FILE  =’ALr,STATUS  =’OLD’) 

OPENCUNTT  =  21, FILE  =’AL2’,STATUS  =’OLD’) 

OPENCUNTT  =  22,FILE  =’AL3’, STATUS  =’OLD’) 

OPENCUNTT  =  23,FILE  = ’OPTAGE’, STATUS  =’OLD’) 

OPENCUNTT  =  24,FILE  = ’LIFETIME’, STATUS  =’OLD’) 

OPENCUNTT  =  25,FILE  = ’MOMENT’, STATUS  =’OLD’) 

OPENCUNTT  =  26,FILE  = ’COUNT’, STATUS  =’OLD’) 

OPENCUNTT  =  27,FILE  =’OPTCOST’, STATUS  =’OLD’) 

CALL  EXCMSC’FILEDEF  50  DISK  FILE  ESTDATl  CDISP  MOD’) 

CALL  EXCMSC’FILEDEF  60  DISK  FILE  ESTDAT2  CDISP  MOD’) 

CALL  EXCMSC’FILEDEF  70  DISK  FILE  ESTDAr3  CDISP  MOD’) 


OPEN(UNIT  =  30, FILE  =’SCALAR2’) 
READ(10,*)OBJ 
IF(OBJ  .GT.  0.05)THEN 
A  =  -1.0 
B  =  0.3 
C  =  0.3 
D  =  0.3 
E  =  -1.0 
F  =  0.3 
G  =  0  3 
H  =  0.3 
I  =  -1.0 
ELSE 

READ(11,*)A 

READ(12,-^)B 

READ(13,*)C 

READ(14,*)D 

READ(15,-^)E 

READ(16,*)F 

READ(17,*)G 

READ(18,*)H 

READ(19,*)I 

ENDIF 

READ(20,*)AL1 
READ(21,*)AL2 
READ(22,’*‘)AL3 
READ(23 ,  ♦)OPTAGE,PI 
READ(24,*)X 

READ(25,*)SECOND, THIRD, AAC, EX 

READ(26,*)COUNT 

READ(27,*)OPTCS 

CLOSE(26) 

PRINT*,’OBJ  =  ’,OBJ 


IF(OPTAGE  .GT.  X)THEN 
Z  =  X 
DI  =  1.0 
ELSE 

Z  =  OPTAGE 
DI  =  0.0 
ENDIF 
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*  OUTPUT 

* 

OPEN(UNIT  =  26, FILE  =’ COUNT’ .STATUS  =’OLD’) 
WRITE(26,  *)COUNT+ 1 
WRITE(50,200)COUNT,A,B,C,D,E,F,G,H,I 
WiaTE(60,300)COUNT,OPTCS,OPTAGE,X,Z,DI,AAC 
WRITE(70, 100)COUNT,AL1  ,AL2,AL3, SECOND, THIRD, OBJ, PI 


*  PREPARE  THE  DATA  FOR  "GAMS”  PROGRAM 
WRITE(30,*)’  SCALAR’ 

WRITE(30,400)’  AA  INITIAL  SOLUTION  OF  A /’, A,’ /’ 

WRITE(30,400)’  BB  INITIAL  SOLUTION  OF  B /’,B,’ /’ 

WRITE(30,400)’  CC  INITIAL  SOLUTION  OF  C  /’,C,’  /’ 

WT?ITE(30,400)’  DD  INITIAL  SOLUTION  OF  D /’,D,’ /’ 

WRITE(30,400)’  EE  INITIAL  SOLUTION  OF  A  /’,E,’  /’ 

WRITE(30,400)’  FF  INITIAL  SOLUTION  OF  B /’,F,’ /’ 
WRITE(30,400)’  GG  INITIAL  SOLUTION  OF  C /’,G,’ /’ 

WRITE(30,400)’  HH  INITIAL  SOLUTION  OF  D  /’,H,’  /’ 

WRITE(30,400)’  H  INITIAL  SOLUTION  OF  D 

WRITE(30,500)’  ALPl  INITIAL  SOLUTION  OF  ALl  /’,AL1,’  /’ 

WRITE(30,500)’  ALP2  INITIAL  SOLUTION  OF  AL2 /’,AL2,’ /’ 

WRITE(30,5(X))’  ALP3  INITIAL  SOLUTION  OF  AL3 /’,AU,’ /’ 

WRITE(30,600)’  MM2  SECOND  STANDARD  MOMENT  /’.SECOND,’/’ 
WRITE(30,600)’  MM3  THIRD  STANDARD  MOMENT /’.THIRD,’/’ 
WRITE(30,700)’  EX  EXPECTED  VALUE  OF  LIFETIME  /’,EX,’  /  ;’ 
100  F0RMAT(1X,I3,3(2X,F5.2),2(2X,F8.5),2X,F4.2,2X,E1 1 .4) 

200  F0RMAT(1X,I3,2X,9(2X,F6.2)) 

300  FORMAT(lX,I3,2X,F8.3,2X,F5.2,2(4X,F7.4),4X,F3. 1 ,4X,F7.3) 

400  FORMAT(A35,F6.2,A2) 

500  FORMAT(A39,F5.2,A2) 

600  FORMAT(A36,F11.8,Al) 

700  FORMAT(A40,F6.3,A4) 

STOP 

END 
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